Wind driven generator bearing state prediction method

Abstract

The invention provides a method for predicting the state of a wind driven generator bearing, and relates to the technical field of electrical equipment fault diagnos.The method comprises the steps that firstly, digital twin models of the wind driven generator bearing under various working conditions are established; the vibration acceleration of each measuring point of the wind driven generator is obtained through electromagnetic-structural field coupling simulation, and a feature database is constructed through obtained signal data; decomposing the signal under the fault condition in the feature database, calculating the fitness of the signal, selecting the optimal signal, and obtaining an optimal modal signal; performing time domain feature extraction on the optimal modal signal; inputting the signal feature vectors measured in different bearing states into the optimized probabilistic neural network, and finally obtaining a model capable of accurately judging the running state of the wind driven generator bearing; vibration acceleration data measured by a detection probe of the wind driven generator in practical application are input into the neural network, the current running state of the wind driven generator bearing can be distinguished, and then monitoring and prediction of the state of the wind driven generator bearing are completed.

Description

Technical Field

[0001] The present invention belongs to the technical field of electrical equipment fault diagnosis, and particularly relates to a method for predicting the bearing state of a wind turbine. Background Art

[0002] During operation, wind turbines generally experience different types of bearing failures due to long-term overload operation and various factors, which can easily lead to a decline in the mechanical performance of the wind turbine and even other structural failures of the wind turbine, affecting the safe and stable operation of the power generation system. At present, the diagnosis of wind turbine bearing faults based on vibration signal analysis has been widely studied. Most of the existing technologies use methods such as wavelet transform and Hilbert-Huang transform to decompose the vibration signal into multiple segments of vibration characteristics, and use machine learning series algorithms to train and diagnose the vibration characteristics and fault types. However, when extracting vibration signal characteristics by the above methods, frequency band aliasing is likely to occur, resulting in the same multiple frequency being repeatedly used as a segmented feature. When performing machine learning diagnosis, the small number of available samples can easily lead to the inability to guarantee the accuracy of the diagnosis algorithm. The diagnosis process lacks a systematic evaluation standard, and the algorithm adaptability is poor. For example, large wind turbines and small wind turbines have different vibration frequency bands from other types, and there are also significant differences in feature distribution and frequency range, resulting in large judgment errors using conventional methods. Existing research shows that the vibration characteristics of the stator and rotor or bearings of motors in severe fault states will exhibit obvious characteristics, but the judgment effect on potential faults such as minor deformations and slight looseness that occur during long-term operation is not good. Moreover, most of the existing technologies are fault diagnosis methods, that is, diagnosing the fault type after the motor fails, and it is impossible to predict the structural abnormalities of the motor before the fault occurs to avoid the occurrence of motor faults. Summary of the Invention

[0003] Object of the Invention

[0004] In order to solve the problems in the existing technology, the present invention provides a method for predicting the hot spot temperature by fusing multi-source temperature information. By using the method of predictive fitting and introducing the BP neural network algorithm optimized by the genetic algorithm to analyze the hot spot temperature rise, the advantages of both the direct measurement method and the indirect measurement method can be combined, and the hot spot temperature of the transformer can be obtained more accurately.

[0005] To achieve the above object, the present invention provides the following technical solutions:

[0006] A method for predicting the bearing state of a wind turbine, comprising the following steps:

[0007] Step 1: Establish a digital twin model of the wind turbine under normal and faulty conditions;

[0008] Step 2: Perform electromagnetic-structural field coupling simulation on the digital twin model of the wind turbine to obtain the electromagnetic force of the stator and the vibration acceleration of the stator;

[0009] Step 3: Analyze the vibration characteristics measured at different measuring points under the first bearing fault condition of the wind turbine based on the digital twin model, and establish a feature database containing the vibration characteristics under the fault condition;

[0010] Step 4: Perform data preprocessing on the data prepared in the database;

[0011] Step 5: Classify the extracted feature values, and judge the state of the first bearing according to the classification results.

[0012] As a further description of the above solution, step 1 includes the following steps:

[0013] Step 1.1: Establish a wind turbine model based on the structural design parameters of the wind turbine generator set and the product drawings;

[0014] Step 1.2: Establish a digital twin model of the wind turbine generator set under three fault conditions of pitting on the outer ring of the first bearing, pitting on the inner ring, and rolling element wear; when establishing the pitting model of the outer ring of the first bearing, change the shape of the outer ring of the first bearing, draw tiny holes on the outer ring to simulate the pitting fault of the outer ring of the first bearing; when establishing the pitting fault of the inner ring of the first bearing, change the shape of the inner ring of the first bearing, draw tiny holes on the inner ring to simulate the pitting fault of the inner ring of the first bearing; when simulating the rolling element wear fault, change the shape of one rolling element of the first bearing, and simulate that the rolling element of the first bearing becomes non-spherical due to wear by drawing a shallow concave shape on one or more rolling elements, so as to simulate the rolling element wear fault of the first bearing; simulate metal fatigue caused by installation error or extrusion by modifying the coaxiality between the first bearing and the bearing housing in the structural long simulation model, add tiny depressions to the outer race to simulate the situation during transportation, modify the friction coefficient between the rolling element and the bearing housing to simulate the situation of bearing sand ingress or other tiny foreign objects, and simulate the situation of rolling element expansion caused by high temperature through the fluid field.

[0015] As a further description of the above solution, step 2 includes the following steps:

[0016] Step 2.1: Perform electromagnetic field finite element simulation according to the model;

[0017] The electromagnetic field simulation is specifically as follows: First, assign the corresponding electromagnetic material properties to each part of the wind turbine. The winding of the physical wind turbine is made of wires. When establishing the digital twin model, a two-dimensional model can be used. In the physical wind turbine, the winding is usually made of wires. When establishing the digital twin model, the winding is regarded as a rectangular planar structure, and the cross-sectional shape of the wind turbine located on the center line of the wind turbine is taken for the rest, as the two-dimensional model of the electromagnetic field simulation. The windings are classified according to the classification of phases, and the number of wires and polarities are set respectively. The setting methods of the type and polarity of the permanent magnet are as follows: When the permanent magnet is wound by multiple wires, the wound type is selected; when the permanent magnet is a solid permanent magnet, the solid type is selected; then set the current expression of the three-phase current; then perform mesh division based on the structural dimension data of the wind turbine generator set, set the total duration and step size of the simulation operation. Finally, cut the tooth part of the stator, separately divide the part with the maximum magnetic density of the tooth part of the stator, and import its electromagnetic force into the structural field part;

[0018] Step 2.2: Using the electromagnetic force per unit area in the air gap obtained from the electromagnetic field simulation calculation of the digital twin model of the wind turbine as the input and loading it on the stator, perform a structural field vibration simulation.

[0019] As a further description of the above solution, the step 3 includes the following steps:

[0020] Step 3.1: Obtain the fault characteristics in the fault state based on the bearing state of the digital twin model of the wind turbine, including: the amplitude of the time-domain vibration signal at multiple measurement points on the stator surface and the vibration transmission situation at different measurement points. When the multi-measurement point vibration signal is converted to the frequency domain, the vibration acceleration amplitude and its change at the fundamental frequency, second harmonic frequency, and their harmonic frequencies;

[0021] Step 3.2: Modify the digital twin model of the wind turbine according to step 1.2 and combine with historical fault data so that it can show the fault state, and construct a simulation fault model for different bearing states of the wind turbine;

[0022] Step 3.3: Perform electromagnetic-structural coupling simulation on the digital twin model of the wind turbine in the normal state and the digital twin model in the fault situation, and extract the vibration accelerations at the measurement point positions p 1, , p 2 , …, p 10 ;

[0023] Step 3.4: Repeat step 3.3 to the set number of times to obtain the vibration accelerations at different measurement point positions as the analysis data samples, and use the vibration acceleration amplitudes including the fundamental frequency and harmonic frequencies in the normal state and the fault state as the frequency domain characteristics of the vibration signal to construct a feature database.

[0024] As a further description of the above solution, step 4 includes the following steps:

[0025] Step 4.1: Randomly initialize the algorithm output parameters, and use this decomposition algorithm to decompose the signal to separate the fault characteristics of the system vibration;

[0026] Step 4.2: Use an optimization algorithm to find the optimal input parameters of the algorithm in step 4.1;

[0027] Step 4.3: Extract time characteristics.

[0028] As a further description of the above solution, step 5 includes the following steps:

[0029] Step 5.1: Integrate the eigenvalue vectors;

[0030] Step 5.2: Use a probabilistic neural network to classify the features;

[0031] Step 5.3: Optimize the input parameters of the probabilistic neural network;

[0032] Step 5.4: Train the probabilistic neural network;

[0033] Step 5.5: Classify the detection data of the wind turbine generator set, and judge the state of the first bearing according to the classification result.

[0034] Advantages and effects of the present invention:

[0035] 1. In step 4 of the present invention, a decomposition algorithm dedicated to mechanical fault diagnosis is adopted, which can extract impact and periodicity, and use filter update and period estimation. Even in the face of complex mechanical faults, it can accurately extract. Compared with the current mainstream algorithms that extract all signal features, steps 4.1 and 4.2 of this algorithm only extract mechanical fault features, so as to intelligently select the features useful for subsequent algorithms. Compared with other current mainstream feature extraction algorithms, it can reduce the calculation pressure of the neural network and improve the accuracy. The decomposition algorithm of the present invention can extract more fault components compared with the traditional decomposition algorithm VMD.

[0036] 2. In steps 1 and 2 of the present invention, a fault prediction method combining digital twin and algorithm is adopted. Even in the face of the problem of fault prediction of equipment with very little effective data, it can accurately obtain the mathematical model for fault prediction, and solve the negative impact caused by insufficient fault data of wind turbines in the field of fault prediction. The current mainstream fault detection methods all rely on manually marked data. Steps 1 and 2 of the present invention can reduce this part of the dependence, have higher universality, and better fill the data types that cannot be obtained by conventional experiments through the finite element method in the digital twin method.

[0037] 3. In step 4.3 of the present invention, a time trend extraction algorithm is used to extract the time trend of the signal. Step 4.3 can retain the time trend as much as possible. Compared with the conventional feature extraction algorithm, step 4.3 of the present invention can better restore the original trend of the data. By extracting the trend of the time series and transforming it into a fluctuation function, the problems of incomplete or over-decomposition in frequency domain decomposition can be avoided.

[0038] 4. Most current mechanical fault prediction methods are based on data obtained from experiments, but experiments will lead to a reduction in economic benefits. In step 3, more types of fault signals obtained based on the digital twin method avoid the dependence on experiments of conventional detection methods and have higher economic value. Through the fluid field analysis of the finite element method, the bearing operation data under high-temperature conditions can be obtained, which also shows that the digital twin method can simulate more bearing operation environments. Therefore, a large amount of environmental data can be used to assist in determining the causes of bearing faults. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 It is a flowchart of a method for predicting the state of a wind turbine bearing according to the present invention;

[0040] Figure 2 It is a vibration equivalent spring-mass model diagram of the wind turbine part of the present invention;

[0041] Figure 3 It is an example electrical engineering drawing of the wind turbine simulation part of the present invention;

[0042] Figure 4 is Figure 3 a schematic structural diagram of the rotor processing housing and the internal rotating shaft in

[0043] Figure 5 is Figure 3 a schematic structural diagram of the rotor in

[0044] Figure 6 is Figure 3 a schematic structural diagram of the stator in

[0045] Figure 7 is Figure 3 a schematic structural diagram of bearing 1 in

[0046] Figure 8 is Figure 3 a schematic structural diagram of bearing 2 in

[0047] Figure 9 is Figure 3 a schematic structural diagram of the gasket in

[0048] Figure 10 is Figure 3 a schematic structural diagram of the motor fixing bracket in

[0049] Figure 11 This is the spring - mass model of the healthy bearing of the present invention;

[0050] Figure 12 This is the spring - mass model of the rolling - element fault bearing of the present invention;

[0051] Figure 13 This is the spring - mass model of the inner - race fault bearing of the present invention;

[0052] Figure 14 This is the spring - mass model of the outer - race fault bearing of the present invention;

[0053] Figure 15 This is the motor vibration transmission diagram of the present invention;

[0054] Figure 16 This is the image of different modal signals obtained by signal decomposition of the present invention;

[0055] Figure 17 This is the effect diagram obtained by training the neural network of the present invention;

[0056] Figure 18 This is the effect diagram for classifying the faults of the wind turbine of the present invention;

[0057] Figure 19 This is the simulation diagram of the example motor of the present invention;

[0058] Figure 20 This is the simulation diagram of the bearing of the example motor of the present invention.

[0059] In the drawings, the list of components represented by each reference numeral is as follows:

[0060] 1 - external rotating shaft, 2 - motor fixing bracket, 3 - rotor processing housing, 4 - rotor, 5 - stator, 6 - internal rotating shaft, 7 - first bearing, 8 - second bearing, 9 - gasket; 10 - outer ring, 11 is the oil film outside the bearing, 12 - rolling element, 13 is the oil film inside the bearing, 14 - inner ring. Detailed implementation manners

[0061] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without making creative efforts belong to the scope of protection of the present invention.

[0062] The present invention proposes a method for predicting the state of a wind turbine bearing. The general flow chart is as shown in Figure 1 and includes the following steps:

[0063] Step 1: Establish digital twin models of the wind turbine under normal and faulty conditions;

[0064] Step 1.1: Based on the structural design parameters of the wind turbine generator set, including the length, width, slot size of the stator 5 punching sheet, the thickness and radius of the punching sheet, the number of punching sheets, the radius of the rotor 4, the rotor thickness, the rotor slot size, the length, width, height and shape of the permanent magnet, the number of permanent magnets, the length and radius of the shaft, the housing size, the motor base size, the motor blade size and shape, the motor bearing model, inner diameter, outer diameter, thickness, rolling element size and number, establish a wind turbine model based on the structural design parameters of each component and the product drawings;

[0065] Such as Figure 3 the external rotating shaft 1, the motor fixing bracket 2, the rotor processing housing 3, the rotor 4, the stator 5, the internal rotating shaft 6, the bearing 1 7, the bearing 2 8, and the gasket 9 in, all need to be modeled according to the structural parameters.

[0066] Step 1.2: Establish digital twin models of the wind turbine generator set under three fault conditions: pitting of the outer ring 10 of the first bearing 7, pitting of the inner ring 14, and wear of the rolling elements; when establishing the pitting model of the outer ring 10 of the first bearing 7, change the shape of the outer ring 10 of the first bearing 7, draw tiny holes on the outer ring 10 to simulate the pitting fault of the outer ring 10 of the first bearing 7; when establishing the pitting fault of the inner ring 14 of the first bearing 7, change the shape of the inner ring 14 of the first bearing 7, draw tiny holes on the inner ring 14 to simulate the pitting fault of the inner ring 14 of the first bearing 7; when simulating the wear fault of the rolling element 12, change the shape of one rolling element 12 of the first bearing 7, by drawing a shallow concave shape on one or more rolling elements, simulate that the shape of the rolling element 12 of the first bearing 7 becomes non-spherical due to wear, so as to simulate the wear fault of the rolling element 12 of the first bearing 7; simulate metal fatigue caused by installation errors or extrusion by modifying the coaxiality between the first bearing 7 and the bearing housing in the structural long simulation model, add tiny depressions to the outer race to simulate the situation during transportation, modify the friction coefficient between the rolling element 12 and the bearing housing to simulate the situation of bearing sand ingress or other tiny foreign objects, and simulate the situation of rolling element expansion caused by high temperature through fluid field simulation. According to Figure 11 the bearing structure in, calculate the outer ring 10, the outer oil film 11, the inner oil film 13, the rolling element 12, and the inner ring 14. If necessary, the oil film can be ignored.

[0067] Step 2: Perform electromagnetic-structural field coupling simulation on the digital twin model of the wind turbine to obtain the electromagnetic force of the stator (4) and the vibration acceleration of the stator (4);

[0068] Step 2.1: Based on the electromagnetic parameters in the product design parameters of the wind turbine, including: rated capacity, rated voltage, rated current, winding winding method, number of permanent magnets, phase current, current density, and the materials selected for the wind turbine structure, conduct electromagnetic field simulation on the digital twin model of the wind turbine;

[0069] Specifically, for the electromagnetic field simulation, first endow each part of the wind turbine structure with corresponding electromagnetic material properties. The winding of the physical wind turbine is made of wires. When establishing the digital twin model, a two-dimensional model can be used; in the physical wind turbine, the winding is usually made of wires. When establishing the digital twin model, the winding is regarded as a rectangular planar structure, and the cross-sectional shape of the wind turbine located on the center line of the wind turbine is taken for the rest, as the two-dimensional model for electromagnetic field simulation; the windings are classified according to the classification of phases, and the number of wires and polarities are set respectively; the setting methods for the type and polarity of permanent magnets are as follows: when the permanent magnet is wound by multiple wires, the wound type is selected; when the permanent magnet is a solid permanent magnet, the solid type is selected; then set the current expression of the three-phase current; then based on the structural dimension data of the wind turbine set, perform mesh division, set the total duration and step size of the simulation operation. Finally, cut the tooth part of the stator 4, separately divide the part with the maximum magnetic density of the tooth part of the stator 4, and import its electromagnetic force into the structural field part;

[0070] Step 2.1.1: When three-phase alternating current is applied to the stator of an asynchronous motor, a series of rotating magnetomotive forces whose factors such as number of times, amplitude, and speed change sinusoidally with time will be generated. The combined magnetomotive force of the stator and rotor windings is:

[0071]

[0072] It can be seen from Equation (1) that the combined magnetomotive force of the stator and rotor windings is mainly divided into three parts, which are respectively

[0073] 2.1.1.1: Main wave combined magnetomotive force

[0074]

[0075] Among them, F 0 is the amplitude of the main wave combined magnetomotive force, ω 1 is the angular frequency of the main wave combined magnetomotive force, is the initial phase angle of the main wave combined magnetomotive force, p is the number of pole pairs, θ is the spatial angle, and t is the time.

[0076] 2.1.1.2: v-th harmonic magnetomotive force of the stator winding

[0077]

[0078] Among them, F vis the amplitude of the synthetic magnetomotive force of the v-th harmonic, ω v is the angular frequency of the synthetic magnetomotive force of the v-th harmonic, is the initial phase angle of the synthetic magnetomotive force of the v-th harmonic.

[0079] 2.1.1.3: μ-th Harmonic Magnetomotive Force of Rotor Winding

[0080]

[0081] Among them, F μ is the amplitude of the synthetic magnetomotive force of the μ-th harmonic, ω v is the angular frequency of the synthetic magnetomotive force of the v-th harmonic, is the initial phase angle of the synthetic magnetomotive force of the μ-th harmonic.

[0082] Since windings are generally three-phase symmetrical windings and the currents fed in are three-phase symmetrical currents, the amplitude of the fundamental wave synthetic magnetomotive force is

[0083]

[0084] Among them, N 1 represents the number of turns of each phase winding, K dp represents the winding coefficient, I m represents the amplitude of the phase current.

[0085] The amplitude of the synthetic harmonic magnetomotive force of the stator three-phase winding is

[0086]

[0087] Among them, K dpv is the winding coefficient of the v-th harmonic, I 1 is the amplitude of the fundamental wave phase current.

[0088]

[0089] Among them, β is the winding intercept coefficient, q is the number of slots per pole per phase, α represents the slot-pole angle, and represents the slot center angle between each phase winding.

[0090] The amplitude of the synthetic harmonic magnetomotive force of the rotor three-phase winding is

[0091]

[0092] Among them, I 2 represents the current of each phase.

[0093] Under ideal conditions, both the stator and rotor 4 of the asynchronous motor have teeth and slots, and the air-gap permeance can be approximately expressed as

[0094]

[0095] Among them, Λ0 represents the fundamental permeance of the air gap, is the k 1 th stator slot harmonic permeance, is the k 2 th rotor slot permeance, is the permeance caused by the interaction between the stator and rotor due to slotting, k 1 represents the harmonic order of the stator slots, k 2 represents the harmonic order of the rotor slots.

[0096] The constant part in the permeance is

[0097]

[0098] where, Λ 0 represents the fundamental permeance of the air gap, μ 0 is the permeability of free space, δ is the thickness of the air gap, K E is the edge effect coefficient of the air gap.

[0099] When the stator has slots and the rotor 4 surface is smooth, the permeance is the stator slot permeance, that is

[0100]

[0101] where, is the amplitude of the k 1 th permeance, which is a constant related to the stator slot structure and material properties, cosk 1 Z 1 θ represents the distribution of the harmonic permeance along the spatial angle θ. This term is the harmonic term, which makes the permeance show periodic variation, Z 1 represents the number of stator slots. This formula can represent the influence of different orders of stator slot harmonics on the permeance.

[0102] When the inner surface of the stator is smooth and the rotor 4 has slots, it is the rotor slot permeance, that is

[0103]

[0104] where, represents the distribution of the harmonic permeance along the spatial angle θ and time t. This term represents the change in permeance caused by the rotation of the rotor 4 and includes the slip s, Z 2 represents the total number of rotor slots, s represents the slip of the rotor 4, which represents the deviation between the rotor 4 and the synchronous speed. ω is the synchronous angular velocity.

[0105] The permeance caused by the interaction between the stator and rotor due to slotting is

[0106]

[0107] where, represents the harmonic magnetic conductance amplitude of the rotor slot, and k is a harmonic coefficient used to represent different harmonic components.

[0108] Introduce Carter coefficient K e , solve and are shown as follows respectively.

[0109]

[0110] Among them, represents the opening coefficient of the stator slot.

[0111]

[0112] Among them, represents the opening coefficient of the rotor slot.

[0113] The calculation formula of the opening coefficient is as follows:

[0114]

[0115] Among them, δ 1 represents the slot opening width, and δ 0 is the tooth pitch.

[0116]

[0117] Among them, b(θ,t) represents the distribution of magnetic induction intensity with respect to the spatial position θ and time t, and v z represents the spatial frequency of the v - th harmonic of the stator winding, and μ Z represents the spatial frequency of the μ - th harmonic of the rotor winding.

[0118] During the operation of the asynchronous motor, there will be a series of electromagnetic force waves in the air gap due to the interaction between the stator and the rotor. According to Maxwell's law, the instantaneous value of the radial electromagnetic force per unit area in the air gap of the motor is

[0119]

[0120] Substitute Equation (17) into Equation (18), and at the same time ignore the force wave components with higher vibration orders but smaller amplitudes and the constant components that have no effect on vibration, we can get

[0121]

[0122] Among them,

[0123] B 1 is the amplitude of the main - wave magnetic flux density, is the phase of the main - wave magnetic flux density, represents the amplitude of the magnetic flux density of the μ - th harmonic of the stator, The amplitude of the magnetic flux density of the v-th harmonic of the rotor 4, ω μ Represents the angular frequency of the μ-th harmonic of the stator, Represents the phase of the μ-th harmonic of the stator, Represents the phase of the 4v-th harmonic of the rotor 4, p r Is the instantaneous value of the radial electromagnetic force.

[0124] Step 2.2: Taking the electromagnetic force loaded on the stator 4 per unit area in the air gap obtained from the electromagnetic field simulation calculation of the digital twin model of the wind turbine as the input, perform a structural field vibration simulation;

[0125] Step 2.2.1: Establish the vibration equation of the stator

[0126] Step 2.2.1.1: Calculate the vibration acceleration generated by the external excitation:

[0127]

[0128] Among them, F 0 Is the amplitude of the electromagnetic force at the spatial position x, ω 0 Is the frequency, i is the imaginary unit, e is the base of the natural logarithm, t represents time, and F(t, x) represents the instantaneous electromagnetic force at the spatial coordinate x and the time coordinate t.

[0129]

[0130] Among them, M is the mass matrix of the stator. The mass matrix describes the mass distribution of the system, and it represents the mass per unit volume. u is the deformation displacement of the stator, C is the damping matrix of the stator. The damping matrix describes the damping effect of the structure on vibration, and it can consider the influence of internal friction of the material and the external damping system. K is the stiffness matrix of the stator.

[0131] Assume that the displacement generated by the external excitation is:

[0132] u ext (t, x) = U ext (x)e iωt (22)

[0133] Among them, u ext (t, x) is the displacement generated by the excitation of the structure in the space-time field, U ext (x) is the distribution of the displacement generated by the excitation of the stator in space, Represents resonance.

[0134] Substituting the spatial part of Equation (21), the spatial distribution of the displacement of the stator generated by the external excitation at the corresponding moment can be obtained:

[0135] U ext(x) = (-Mω 0 2 + iω 0 C + K) -1 F 0 (x) (23)

[0136] Among them, since the wind turbine is a pure rigid body structure, the above formula can be simplified to:

[0137] U ext (x) = (-Mω 0 2 + K) -1 F 0 (x) (24)

[0138] Taking the second-order derivative of Equation (24) and taking the absolute value can obtain the vibration acceleration amplitude of the stator generated by the external excitation at the corresponding moment:

[0139]

[0140] Among them, represents the absolute value of the vibration acceleration of the stator under the action of the external excitation, and |U ext (x)| represents the deformation displacement of the stator under the action of the external excitation.

[0141] Step 2.2.1.2: Calculate the vibration acceleration under the natural frequency vibration

[0142] Determine the natural frequency of the stator through eigenvalue analysis. In free vibration, considering the case of no external load, the dynamic equilibrium equation is:

[0143]

[0144] To find the natural frequency and the corresponding mode, assume the solution is:

[0145]

[0146] where φ is the mode shape and ω n is the natural frequency.

[0147] Substitute the solution of Equation (26) into Equation (25):

[0148] The first term:

[0149]

[0150] The second term:

[0151]

[0152] After substitution, Equation (28), the dynamic equilibrium equation becomes:

[0153]

[0154] Further simplifying Equation (29) gives:

[0155] (-Mω n 2 +iω n C + K)φ = 0 (30) Finally, the eigenvalue simplification problem is obtained:

[0156] det(-Mω n 2 +iω n C + K) = 0 (31)

[0157] In a wind turbine, since all components are rigid bodies, the damping matrix C is a zero matrix, and we get:

[0158]

[0159] where ω n is the natural frequency to be solved.

[0160] Solve for the mode shape corresponding to the natural frequency:

[0161]

[0162] where φ n is the mode shape corresponding to the natural frequency ω n

[0163] The relationship between displacement and mode shape is as follows:

[0164] v nat (t, x) = φ(x)q(t) (34)

[0165] where q(t) is the time response, and u n (t, x) is the displacement of the natural vibration, which is usually in the form of a harmonic vibration:

[0166]

[0167] where A is the amplitude.

[0168] Taking the second derivative of Equation (35) gives:

[0169]

[0170] Substituting Equation (36) into Equation (25) gives:

[0171]

[0172] ​Substituting Equation (36) into Equation (37) gives:

[0173]

[0174] Taking the real part of Equation (38) gives the vibration acceleration generated by the natural frequency vibration:

[0175]

[0176] Step 2.2.1.3: Calculate the vibration acceleration of the stator of the wind turbine under the action of an external electromagnetic force excitation

[0177] The calculation formula for the vibration acceleration generated by a rigid body under the action of an external excitation is as follows:

[0178]

[0179] where represents the vibration acceleration generated by the external excitation, represents the vibration acceleration generated by the natural frequency of the rigid body. is the total vibration acceleration generated when the stator of the wind turbine is subjected to an external excitation. In Steps 1 and 2 of the present invention, a fault prediction method combining digital twin and algorithms is adopted. Even for the fault prediction problem of equipment with very little valid data, the mathematical model for fault prediction can be accurately obtained, solving the negative impact caused by insufficient fault data of wind turbines in the field of fault prediction. The current mainstream fault detection methods all rely on manually marked data. Steps 1 and 2 of the present invention can reduce this part of the dependence and have higher universality. As shown in the appendix Figure 17 , the finite element method in the digital twin method is used to better fill in the data types that cannot be obtained by conventional experiments.

[0180] Step 3: Analyze the vibration characteristics measured at different measurement points under the bearing fault condition of the wind turbine based on the digital twin model, and establish a feature database containing the vibration characteristics under the fault condition;

[0181] Step 3.1: Obtain the fault characteristics under the fault state based on the bearing state of the wind turbine digital twin model, including: the amplitude of the time-domain vibration signal at multiple measurement points on the surface of the stator 4 and the vibration transmission situation at different measurement points, and the amplitude and change situation of the vibration acceleration at the fundamental frequency, second harmonic frequency and their multiple frequencies when the multi-measurement point vibration signal is converted to the frequency domain;

[0182] Step 3.2: Modify the wind turbine digital twin model according to Step 1.2 and in combination with historical fault data so that it can exhibit the fault state, and construct a simulation fault model for different bearing states of the wind turbine;

[0183] Step 3.3: Conduct electromagnetic-structural coupling simulation on the digital twin model of the wind turbine in normal state and the digital twin model under fault conditions, and extract the vibration accelerations at the measuring point positions p 1, , p 2 , …, p 10 under normal operation and fault conditions.

[0184] Step 3.4: Repeat Step 3.3 for a set number of times to obtain the vibration accelerations at different measuring point positions as the analysis data samples. Take the vibration acceleration amplitudes including the vibration fundamental frequency and harmonic frequencies under normal and fault conditions as the frequency domain characteristics of the vibration signals and construct a feature database. Currently, most mechanical fault prediction methods are based on data obtained from experiments, but experiments can lead to reduced economic benefits. In Step 3, more types of fault signals obtained based on the digital twin method avoid the dependence on experiments of conventional detection methods and have higher economic value. As Figure 19 shown, through the fluid field analysis of the finite element method, the bearing operation data under high-temperature conditions can be obtained, which also shows that the digital twin method can simulate more bearing operation environments. Therefore, a large amount of environmental data can be used to assist in determining the causes of bearing faults.

[0185] Step 4: Perform data preprocessing on the data prepared in the database. The data preprocessing process is as follows:

[0186] Step 4.1: Randomly initialize the algorithm output parameters, and use this decomposition algorithm to decompose the signal to separate the fault characteristics of the system vibration;

[0187] Step 4.1.1: Load the original signal, input the parameters filter length L and number of modes K (the values of L and K are confirmed through Step 4.2), and mark the input signal of length N as X(N), expressed as:

[0188] X p (N) = {x 1 , x 2 , …, x n}, p = 1, 2, 3, …, 10; N = 1, 2, 3, …, N (41)

[0189] where p represents the index of the detection probe, and x n represents the vibration acceleration amplitude at the corresponding moment, and N is the sequence length. Step 4.1.2: Initialize the FIR filter bank through the Hanning window, use K filters (generally recommended 5 - 10), set the iteration coefficient i = 1, and start the iteration.

[0190] The calculation method of the FIR filter coefficients is as follows:

[0191] 4.1.2.1: First, divide the frequency band of the original signal into K segments to obtain the cut-off frequencies:

[0192]

[0193] where: f s is the sampling frequency of the original signal, f l is the lower cut-off frequency, and f u is the upper cut-off frequency.

[0194] 4.1.2.2: Initialize the FIR filter

[0195] The coefficients of the ideal high-pass FIR filter are calculated as follows:

[0196]

[0197] The coefficients of the ideal low-pass FIR filter are calculated as follows:

[0198]

[0199] The coefficients of the ideal band-pass FIR filter are:

[0200] f(l) = f s (l) - f l (l) (45)

[0201] After initialization through the Hanning window:

[0202] f(l) = (f s (l) - f l (l))iω L (46)

[0203] where ω L represents the Hanning window, and its expression is as follows:

[0204]

[0205] 4.1.2.3: Start the iteration. To obtain the filter coefficients that can maximize the relevant kurtosis, convert the optimization problem into a constrained problem:

[0206] Adopt a:

[0207]

[0208] where u k represents the k-th decomposition mode, f k is the coefficient of the k-th filter, its length is L, T s is the sampling period, M is the shift order, CK MDenote the relevant kurtosis corresponding to the decomposition mode, and use it as the objective function. The sampling period is T s is the period of this kind of fault, and the measurement method of this value will be explained in step 4.1.4

[0209] 4.1.2.4: Use an iterative eigenvalue decomposition method to solve the constrained problem in equation (48)

[0210] u k = Xf k (49)

[0211] where:

[0212]

[0213] 4.1.2.5: Transform equation (48) to obtain the optimization problem

[0214] Substitute equation (49) into equation (48) to get:

[0215]

[0216] where, denotes the conjugate transpose operation on u k and W M is a control weighting matrix, and its content is as follows:

[0217]

[0218] where:

[0219]

[0220] Substitute equation (49) into equation (48) to obtain:

[0221]

[0222] where, R XWX is the weighted correlation matrix, and R XX is the correlation matrix. R XX The calculation method is as follows:

[0223]

[0224] 4.1.2.6: Equivalent the optimization problem in equation (47) to the problem of finding eigenvectors

[0225] Maximizing the filter coefficients by solving equation (48) can be transformed into generally solving the following formula to obtain the eigenvector corresponding to the maximum eigenvalue λ:

[0226] R XWX f k = R XXf k λ (56)

[0227] Therefore, the coefficients of the k-th filter will be updated by Equation (56), and iterated continuously to approximate the set target, that is, the filtered signal with the maximum correlation kurtosis.

[0228] Step 4.1.3: Calculate the signal of the k-th mode in the i-th iteration through the following formula

[0229]

[0230] where represents the signal of the k-th mode in the i-th iteration, The coefficient vector of the filter of the signal of the k-th mode in the i-th iteration, * represents the convolution operation, and its calculation expression is as follows:

[0231]

[0232] where represents the n-th signal of the k-th mode in the i-th iteration, represents the l-th filter coefficient of the k-th mode in the i-th iteration, and (n - l) is the shift order M.

[0233] Step 4.1.4: Update the filter coefficients

[0234] Use the original signal X(N), the decomposed modal signal and use the estimated period T k i as The point where the autocorrelation spectrum reaches a local maximum after passing through the zero crossing Update the filter coefficients and complete one iteration, i = i + 1.

[0235] 4.1.4.1: The autocorrelation spectrum of the signal will produce a local maximum at the period position. If R x (τ) represents the autocorrelation function of the signal X(N), and the expression of R x (τ) with respect to the lag τ can be defined as

[0236]

[0237] 4.1.4.2: According to this autocorrelation spectrum function, find the point τ where the local maximum is reached after passing through the zero crossing of the function 1 as the input period, that is:

[0238] T s = τ 1 (60)

[0239] Step 4.1.5: Determine whether the pre-iteration count has been reached. If the condition is met, proceed to the next step; otherwise, return to Step 4.1.3. Step 4.1.6: Calculate the correlation coefficient r between every two modes, construct a K×K matrix, find the two modes with the largest correlation coefficient, and use T s Calculate their correlation kurtosis CK M (u k ), then select the mode with the smaller correlation kurtosis and discard it, and set K = K - 1.

[0240] The calculation expression of the correlation coefficient r is as follows:

[0241]

[0242] where represents the nth signal of mode k 1 , represents the nth signal of mode k 2 , represents the mean of all signals of mode k 1 , represents the mean of all signals of mode k 2 .

[0243] The expression for calculating the mean of a mode is as follows:

[0244]

[0245] Step 4.1.7: Determine whether the number of modes K can reach the preset value. If it cannot be reached, return to Step 4.1.3; otherwise, proceed to the next step.

[0246] Step 4.1.8: Use the remaining modes as the result of the final mode decomposition. The mode result can be expressed as:

[0247]

[0248] where Mode represents the modal component.

[0249] Step 4.2: Use an optimization algorithm to find the optimal input parameters of the algorithm in Step 4.1;

[0250] Step 4.2.1: Set the numerical range of K to (3, 8) and the numerical range of L to (5, 10). Calculate the fitness function for the decomposed modes obtained from each combination of individual input parameters. The calculation method of the fitness function is as follows:

[0251] Step 4.2.1.1: Construct the required feature vector group, including the following steps:

[0252] 4.2.1.1.1: Construct the coarse-grained form of each component scale:

[0253] The coarse-grained calculation process for any Mode is as follows:

[0254]

[0255] Among them, τ is the time scale, usually taken as (5, 10), and it needs to be divisible by N, and the value must be determined according to the empirical method; j is the coarse-grained sequence index, and n is the measured data point index.

[0256] 4.2.1.1.2: Calculate the energy of all coarse-grained time series of Mode k at the corresponding scale factor τ

[0257]

[0258] 4.2.1.1.3: Convert all of Mode k to vector form:

[0259]

[0260] 4.2.1.1.4: The eigenvectors E of all Modes of X(N) can be represented as a vector group: k

[0261] E = [E 1 , E 2 , …, E k , …, E K T , k = 1, 2, …, K (67)

[0262] Step 4.2.1.2: Calculate the sample entropy of E

[0263] The sample entropy calculation process is as follows:

[0264] 4.2.1.2.1: Select the embedding dimension α' and the tolerance β'. α' must be determined according to the empirical method and is usually taken as (2, 10). The calculation expression for the tolerance β' is as follows:

[0265] β' = c × σ (68)

[0266] Among them, c is the proportionality coefficient, usually taken as (0.1, 0.2). σ is the standard deviation of the sequence, and its calculation expression is as follows:

[0267]

[0268] 4.2.1.2.2: Generate all subsequences of length α'. For the sequence E k ​​​​, the j-th subsequence is:

[0269]

[0270] 4.2.1.2.3: Calculate the matching pair subsequences:

[0271] For each pair of subsequences and calculate the distance between them. The calculation expression is as follows:

[0272]

[0273] 4.2.1.2.4: Calculate the number of matching pairs:

[0274] For each pair of subsequences and If then count it as a matching pair and make the following definition: O α； O(β') represents the number of similar pair subsequences among all subsequences of length α, O α'+1 O(β') represents the number of similar pair subsequences among all subsequences of length α'+1.

[0275] 4.2.1.2.5: Calculate the sample entropy

[0276]

[0277] Among them, represents the sample entropy

[0278] Step 4.2.1.3: Obtain the fitness function value:

[0279]

[0280] Step 4.2.2: Calculate the comprehensive fitness function value:

[0281] Step 4.2.2.1: Construct a multi-measurement point fitness function matrix for each combination of K and L:

[0282]

[0283] Among them, δ 1 K represents the K-th mode decomposed by the first measurement point.

[0284] Step 4.2.2.2: Nondimensionalization processing:

[0285]

[0286] Among them, δ′ ij is the data after nondimensionalization processing, δij represents the value at the i-th row and j-th column, max represents finding the maximum value in a set of vectors, and min represents finding the minimum value in a set of vectors.

[0287] Step 4.2.2.3: Calculate the index variability:

[0288]

[0289] is the average value of the j-th column.

[0290] Index variability S j is:

[0291]

[0292] Step 4.2.2.4: Calculate the index conflict R j :

[0293]

[0294] where r ij is calculated as follows:

[0295]

[0296] Step 4.2.2.5: Calculate the information content C of the j-th column j :

[0297] C j = S j × R j (80)

[0298] Step 4.2.2.6: Assign the objective weight W to the j-th column j :

[0299]

[0300] Step 4.2.2.7: Complete the calculation of the comprehensive fitness function:

[0301]

[0302] where ω' is the value of the comprehensive fitness function at different measurement points.

[0303] Step 4.2.2.8: Complete the screening of the comprehensive fitness function:

[0304] Ω = min{ω' 1 , ω' 2 , …, ω' 10} (83)

[0305] Among them, Ω is the fitness function value under the current combination of K and L.

[0306] Step 4.2.3: Select K and L corresponding to Ω under the combination of the minimum parameters K and L as input parameters.

[0307] Step 4.3: Extract time features, including the following steps:

[0308] Step 4.3.1: Assume the decomposed modes are as follows:

[0309]

[0310] Step 4.3.2: Calculate the mean of the modal components:

[0311]

[0312] Step 4.3.3: Construct a new sequence, including the following steps:

[0313] 4.3.3.1: Construct a new sequence θ k (n):

[0314]

[0315] 4.3.3.2: Group θ k (n):

[0316] n s = [n / s] (87)

[0317] Among them, [n / s] represents the integer operation.

[0318] Step 4.3.4: Fit the fluctuation trend of each sub-interval:

[0319] Use polynomial fitting to fit the trend of each sub-interval to obtain the fluctuation trend y f (n):

[0320]

[0321] Among them, a r' is the coefficient of the r'-order fitting polynomial.

[0322] Step 4.3.5: Eliminate the fluctuation trend of each sub-interval:

[0323] △y f (n) = y(n) - y f (n) (89)

[0324] Step 4.3.6: Calculate the mean square fluctuation value within each sub-interval:

[0325]

[0326] Among them, τ' is the sub - interval index and s' is the number of intervals.

[0327] Step 4.3.7: Calculate the second - order wave function F q (s'):

[0328]

[0329] Step 4.3.8: Modify the sub - interval interval length in Step 4.3.3, and repeat Steps 4.3.3 to 4.3.6 to obtain the fluctuation function F q (s'), and use the least - squares method to fit the data points through the least - squares function to obtain a linear function.

[0330] log 10 F q (s') = H α” log 10 s' + log 10 A α” (92)

[0331] Among them, H α” is the extracted eigenvalue. In Step 4 of the present invention, a decomposition algorithm dedicated to mechanical fault diagnosis is adopted, which can extract impact and periodicity, and utilize filtering update and period estimation. Even in the face of complex mechanical faults, accurate extraction can be achieved. Compared with the current mainstream algorithms that extract all signal features, Steps 4.1 and 4.2 of this algorithm only extract mechanical fault features, thus intelligently selecting features useful for subsequent algorithms. Compared with other current mainstream feature extraction algorithms, it can reduce the computational pressure of the neural network and improve the accuracy. Refer to the appendix Figure 16 It can be seen that the decomposition algorithm of the present invention can extract more fault components compared with the traditional decomposition algorithm VMD; meanwhile, in Step 4.3, a time - trend extraction algorithm is adopted to extract the time trend of the signal. Step 4.3 can retain the time trend as much as possible. Compared with the conventional feature extraction algorithms, Step 4.3 of the present invention can more restore the original trend of the data. Refer to the appendix Figure 18 , by extracting the trend of the time series and transforming it into a fluctuation function, the problems of incomplete frequency - domain decomposition or over - decomposition can be avoided.

[0332] Step 5: Classify the extracted eigenvalues and judge the bearing state according to the classification results, which specifically include the following steps: Step 5.1: Integrate the eigenvalue vectors:

[0333] The features extracted from each signal are in the following form:

[0334] Feature = {H 1 ,H2 ,…,H K} (93)

[0335] Among them, K is the number of modes

[0336] Step 5.2: Classify the features using a probabilistic neural network, including the following steps:

[0337] 5.2.1: Each pattern unit forms the dot product of the input pattern vector Feature and the weight vector W', Z mn Then perform a non-linear operation before outputting its activation level to the summing unit mn = FeatureiW' mn After normalizing Feature and W' σ' is the scaling factor. The PNN decision boundary is associated with the smoothing factor. After normalization, the non-linear operation of Z mn is converted to: mn

[0338]

[0339] where W mn represents the nth training vector of the mth class, and T represents the transpose.

[0340] 5.2.2: Calculate by the neurons in the sample layer:

[0341]

[0342] where p is the dimension of the feature vector, m is the class number, and n is the pattern number.

[0343] 5.2.3: Calculate by the neurons in the summing layer. The clustering pattern Feature is classified as C m (m = 1, 2, 3... N). Calculate the output of the comprehensive calculation neuron:

[0344]

[0345] where N m is the number of training patterns.

[0346] 5.2.4: The decision rule of Bayes is d(x) = C m .

[0347] If:

[0348] l m p(C m )p m (Feature|C m ) ≥ l k p(C k )p k ​(Feature|C k ), k≠m (97)

[0349] where l m is the loss of decision-making error for each category C m , P(C m ) is the probability that the vector conforms to the category, and p m (Feature|C m ) is the conditional probability density of Feature for the class. Classify the cluster pattern Feature according to the Bayesian decision rule:

[0350]

[0351] where is the evaluation classification of the feature Feature, and N is the total number of classes of training samples.

[0352] Step 5.3: Optimize the input parameters of the probabilistic neural network

[0353] To find the optimal value of the smoothing factor σ of the probabilistic neural network, the following method is used for calculation:

[0354] 5.3.1: Set the input parameters of the optimization algorithm: The optimization algorithm needs to input the following three parameters: MaxIter, SearchAgents, and (lb, ub). Generally speaking, MaxIter takes (20, 40), and SearchAgents takes (10, 30). (lb, ub) are the upper and lower bounds of the smoothing factor σ' to be optimized, and here (0.001, 1000) is taken. The fitness function is the classification accuracy of the probabilistic neural network.

[0355] 5.3.2: Calculate the local optimal solution:

[0356] Randomly extract SearchAgents parameters from (lb, ub):

[0357] {σ' 1 , σ' 2 , …, σ' SearchAgents} (99)

[0358] Use them as input parameters and input them into the probabilistic neural network to obtain SearchAgents groups of classification accuracy results:

[0359] {ar” 1 , ar” 2 , …, ar” SearchAgents} (100)

[0360] where ar” 1Represents the accuracy value of the first classification result.

[0361] Find the values of σ corresponding to the largest three of these values as the proxies for α", β", and δ".

[0362] 5.3.3: Calculation expression for the distance between the search agent SearchAgents and the optimal solution As follows:

[0363]

[0364] Among them, is the coefficient vector, is the position of the optimal solution after t iterations, is the position of the search agent after t iterations. The calculation expression of

[0365]

[0366] Among them, is a vector randomly generated between [0, 1].

[0367] 5.3.4: The expression for the position of the updated search agent is as follows:

[0368]

[0369] Among them, represents the position of the search agent after one iteration update, is another coefficient vector, and its calculation expression is as follows:

[0370]

[0371] Among them, is a vector randomly generated between [0, 1]. is the convergence factor, which linearly decreases from 2 to 0 as the iteration progresses. 5.3.5: The mathematical model for the advancing direction of the search agent is as follows:

[0372]

[0373] Among them, respectively represent the distances between the α", β", and δ" proxies and other agents, is the coefficient vector.

[0374]

[0375] Among them, is the advancing direction and step size of the search agent each time, respectively represent the current positions of the α", β", and δ" proxies, A1 , A 2 , A 3 are coefficient vectors respectively.

[0376] 5.3.6: The mathematical model of the finally found optimal solution is as follows:

[0377]

[0378] where is the updated position. Through continuous iteration, the optimal solution is finally found.

[0379] Step 5.4: Train the probability neural network

[0380] Divide the labeled samples (usually the data obtained by simulation) into a training set and a test set (generally allocated in a ratio of 0.8:0.2). Extract features from them through Step 4, and then input them into the probability neural network optimized by Step 5.3 for training.

[0381] Step 5.5: Classify the detection data of the wind turbine generator set, and judge the state of the first bearing 7 according to the classification result. Specifically:

[0382] Input the vibration acceleration data measured by the detection probe on the wind turbine generator during on-site inspection into Steps 4 to 5 above to obtain the state of the bearings of the current wind turbine generator. This method can classify the operating states of the wind turbine generator bearings into four categories: healthy state, inner ring 14 fault, outer ring 10 fault, and rolling element 12 fault.

[0383] Figure 16 is the different modal signal images obtained by signal decomposition in Step 4 of the present invention;

[0384] Figure 17 is the effect diagram obtained by training the neural network in Step 5.4 of the present invention, with an error of 0, and all samples are classified correctly;

[0385] Figure 18 is the effect diagram of classifying the wind turbine generator faults in Step 5.5 of the present invention. All test groups are classified correctly, and the classification accuracy rate is 100%;

[0386] In the actual use of the present invention, based on the actual operating conditions of the wind turbine, the vibration acceleration signals near the rotor and the rotating shaft of the wind turbine are collected through sensors and acquisition devices. The vibration signals are decomposed to obtain multiple modal components. The fitness function of the modal components is calculated, and the decomposition method with the lowest fitness function value is selected. Weight distribution is performed on the multiple modal components obtained by this decomposition method for multiple measurement points, and time features are extracted for each mode. The time features of each mode are combined to form the comprehensive time feature of the signal. Finally, classification is performed through an optimized probabilistic neural network to realize the prediction of the bearing state of the wind turbine generator set.

[0387] As described above, it is only an embodiment of the present application and does not impose any form of limitation on the present application. Although the present application is disclosed as a preferred embodiment, it is not intended to limit the present application. Any person skilled in the art, without departing from the scope of the technical solution of the present application, makes some changes or modifications using the technical content disclosed above, which are equivalent to equivalent implementation cases and all fall within the scope of the technical solution.

Claims

1. A method for predicting the bearing status of a wind turbine generator set, characterized in that: The following steps are involved: Step 1: Establish digital twin models of wind turbines in normal and fault conditions; Step 2: Perform electromagnetic-structural field coupling simulation on the digital twin model of the wind turbine to obtain the electromagnetic force of the stator (4) and the vibration acceleration of the stator (4); Step 3: Analyze the vibration characteristics measured at different measuring points under the fault condition of the first bearing (7) of the wind turbine based on the digital twin model, and establish a characteristic database containing the vibration characteristics under the fault condition; Step 4: Preprocess the data prepared in the database; Step 5: Classify the extracted characteristic values ​​and determine the state of the first bearing (7) according to the classification result.

2. The method for predicting the bearing state of a wind turbine generator set according to claim 1, characterized in that: The step 1 comprises the following steps: Step 1.1: Establish a wind turbine model based on the wind turbine structure design parameters and product drawings; Step 1.2: Establish a digital twin model of the wind turbine generator set under three fault conditions: pitting corrosion of the outer ring (10) of the first bearing (7), pitting corrosion of the inner ring (14), and rolling element wear; when establishing the pitting corrosion model of the outer ring (10) of the first bearing (7), change the shape of the outer ring (10) of the first bearing (7), draw tiny holes on the outer ring (10), and simulate the pitting corrosion fault of the outer ring (10) of the first bearing (7); when establishing the pitting corrosion fault of the inner ring (14) of the first bearing (7), change the shape of the inner ring (14) of the first bearing (7), draw tiny holes on the inner ring (14), and simulate the pitting corrosion fault of the inner ring (14) of the first bearing (7); When simulating the wear failure of the rolling element (12), the shape of a rolling element (12) of the first bearing (7) is changed. By drawing a shallow concave shape on one or more rolling elements, it is simulated that the rolling element (12) of the first bearing (7) of the shaft becomes non-spherical due to wear, thereby simulating the wear failure of the rolling element (12) of the first bearing (7); by modifying the coaxiality of the first bearing (7) and the bearing chamber in the structural length simulation model, the metal fatigue caused by installation error or extrusion is simulated, and the outer rolling ring is increased with a small depression to simulate the situation during transportation, and the friction coefficient between the rolling element (12) and the bearing chamber is modified to simulate the situation of sand or other small foreign matter entering the bearing, and the rolling element expansion problem caused by high temperature is simulated by the fluid field.

3. The method for predicting the bearing state of a wind turbine generator set according to claim 1, characterized in that: The step 2 comprises the following steps: Step 2.1: Perform electromagnetic field finite element simulation based on the model; The electromagnetic field simulation is specifically as follows: first, the electromagnetic material properties corresponding to the structure of each part of the wind turbine are assigned. The winding of the physical wind turbine is wound by a wire. When establishing the digital twin model, a two-dimensional model can be used. In the physical wind turbine, the winding is usually wound by a wire. When establishing the digital twin model, the winding is regarded as a rectangular plane structure, and the rest of the parts are the cross-sectional shape of the wind turbine located at the center line of the wind turbine as the two-dimensional model of the electromagnetic field simulation; the winding is classified according to the phase classification, and the number of wires and polarity are set respectively; the method for setting the type and polarity of the permanent magnet is as follows: when the permanent magnet is wound by multiple wires, the winding type is selected; when the permanent magnet is a permanent magnet solid, the solid type is selected; then the current expression of the three-phase current is set; then the grid is divided based on the structural size data of the wind turbine generator set, the total duration and step length of the simulation run are set, and finally, the teeth of the stator (4) are cut, and the part with the largest magnetic density of the teeth of the stator (4) is separated separately, and its electromagnetic force is introduced into the structural field part; Step 2.2: The electromagnetic force per unit area in the air gap obtained by the electromagnetic field simulation calculation based on the digital twin model of the wind turbine and loaded on the stator (4) is used as input to perform structural field vibration simulation.

4. The method for predicting the bearing state of a wind turbine generator set according to claim 1, characterized in that: The step 3 comprises the following steps: Step 3.1: Obtain fault characteristics under fault conditions based on the bearing status of the digital twin model of the wind turbine, including: the amplitude of the time-domain vibration signal at multiple measuring points on the surface of the stator (4) and the vibration transmission conditions at different measuring points, and the amplitude and change of the vibration acceleration at one times the fundamental frequency, two times the fundamental frequency and their multiples when the vibration signals at multiple measuring points are converted to the frequency domain; Step 3.2: According to step 1.2 and combined with historical fault data, modify the wind turbine digital twin model so that it can show the fault state and build a simulation fault model of different wind turbine bearing states; Step 3.3: Perform electromagnetic-structural coupling simulation on the digital twin model of the wind turbine in normal state and the digital twin model in fault state to extract the measurement point positions p under normal operation and different fault states. 1, ,p2,…,p 10 Vibration acceleration; Step 3.4: Repeat step 3.3 to the set number of times to obtain the vibration acceleration at different measuring points as analysis data samples, and use the vibration acceleration amplitude including the vibration fundamental frequency and multiple frequency under normal and fault conditions as the frequency domain characteristics of the vibration signal and build a feature database.

5. The method for predicting the bearing status of a wind turbine generator set according to claim 1, characterized in that: The step 4 comprises the following steps: Step 4.1: Randomly initialize the algorithm output parameters, and use the decomposition algorithm to decompose the signal to separate the fault characteristics of system vibration; Step 4.2: Find the optimal input parameters of the algorithm in step 4.1 through the optimization algorithm; Step 4.3: Extract temporal features.

6. The method for predicting the bearing status of a wind turbine generator set according to claim 1, characterized in that: The step 5 comprises the following steps: Step 5.1: Integrate the eigenvalue vector; Step 5.2: Use probabilistic neural network to classify features; Step 5.3: Optimize the probabilistic neural network input parameters; Step 5.4: Train the probabilistic neural network; Step 5.5: Classify the wind turbine generator set detection data, and determine the state of the first bearing (7) according to the classification result.

Images