A wind farm wind turbine intelligent monitoring method

Abstract

The application relates to the technical field of intelligent monitoring, and discloses a wind farm fan intelligent monitoring method, which comprises the following steps: under the same sampling clock reference and synchronous time stamp, synchronously reading a generator time domain active power sequence and a time domain rotating speed sequence, and obtaining a space angle domain equivalent torque sequence with a wind wheel mechanical azimuth as an independent variable through integration and resampling; extracting a single-circle average reference torque and combining a cosine analytical function to synthesize a three-fold frequency reference aerodynamic excitation sequence; performing continuous cross-correlation integral on the two to analyze a transmission chain torsional hysteresis phase angle; generating a dynamic equivalent torsional stiffness characteristic value according to the ratio of the integral average of the current cycle objective work torque and the transmission chain torsional hysteresis phase angle, solving a smoothing bandwidth according to the standard deviation to construct a state probability density function and generate system state information entropy; and finally deducing a health state deviation dispersion through cross-logarithmic integral transformation, and converting the health state deviation dispersion into a failure early warning probability sequence in combination with a historical degradation variation factor.

Description

Technical Field

[0001] This invention relates to the field of intelligent monitoring technology, and more specifically, to an intelligent monitoring method for wind turbines in a wind farm. Background Technology

[0002] Megawatt-class wind turbines are typically equipped with flexible composite blades ranging from hundreds of meters in length. The rear end of the rotor is connected to a multi-body torsional flexible transmission chain, including a low-speed main shaft, a multi-stage gearbox, and a high-speed shaft. In natural operating environments, the atmospheric boundary layer of the wind farm follows a non-uniform natural wind shear distribution in vertical space, and the airflow passing over the lower cylindrical tower generates a tower shadow effect. Due to these natural conditions, the turbine components inevitably experience periodic alternating aerodynamic loads during the rotational sweep process.

[0003] Under prolonged full-load operation, the internal materials of the transmission chain are prone to microscopic fatigue degradation, leading to an irreversible decrease in its torsional equivalent stiffness. Current monitoring and diagnostic systems, when processing massive amounts of asynchronous operating data, typically treat data within a specified time window as a rigid, synchronous evolution state, failing to effectively address the spectral dispersion effect caused by variable speed. Due to the highly random nature of turbulent wind speeds in nature, the power and vibration amplitude collected by the generator set are constantly oscillating within a wide frequency range. This complex background, composed of macroscopic aerodynamic noise, easily masks weak variation signals caused by the degradation of the machine's stiffness.

[0004] Therefore, existing conventional retrieval and monitoring algorithms are prone to misinterpreting drastic changes in external environmental conditions as abnormalities in internal mechanical components, resulting in a high false alarm rate. On the other hand, when the machine body actually shows signs of microscopic fatigue fracture, the warning indicators are easily drowned out by macroscopic environmental noise, leading to missed reports and making it difficult to trace and accurately locate early hidden dangers before irreversible damage to components. Summary of the Invention

[0005] This invention provides an intelligent monitoring method for wind turbines in wind farms, which solves the technical problems mentioned in the background art.

[0006] This invention provides a method for intelligent monitoring of wind turbines in a wind farm, applied to an electronic digital data processing system including a central processing unit, comprising:

[0007] Under the same sampling clock reference and synchronization timestamp, the generator time-domain active power sequence and time-domain speed sequence are read synchronously. After continuous integration and spline interpolation resampling, the spatial angular domain equivalent torque sequence with the wind turbine mechanical azimuth angle as the independent variable is obtained.

[0008] Within the azimuth domain of a single wind turbine rotation, with the mechanical angles from 0 to 2π corresponding to one revolution of the wind turbine as the range and the azimuth angle corresponding to the Z pulse of the main shaft encoder as the starting point, the integral expectation of the equivalent torque sequence in the spatial azimuth domain is extracted as the average reference torque for a single revolution, and combined with the third harmonic cosine analytical function to synthesize the third harmonic reference aerodynamic excitation sequence.

[0009] By performing cross-correlation integration on the third harmonic reference aerodynamic excitation sequence and the spatial angular domain equivalent torque sequence, the torsional hysteresis phase angle of the transmission chain is analyzed.

[0010] The ratio of the average integral value of the objective working torque of the current cycle to the torsional hysteresis phase angle of the transmission chain is used as the current dynamic equivalent torsional stiffness characteristic value. The current dynamic equivalent torsional stiffness characteristic value is used to characterize the average torque carrying capacity corresponding to the unit torsional hysteresis phase angle of the transmission chain and is collected into a stiffness characteristic sample set.

[0011] The optimal smoothing bandwidth is calculated based on the discrete standard deviation of the stiffness feature sample set, and the state probability density function is constructed and integrated to generate the system state information entropy.

[0012] Retrieve the initial health steady-state baseline probability distribution, perform a cross-logarithmic integral transformation on the state probability density function, deduce the current health state deviation divergence and update it to the historical deviation divergence set;

[0013] The squared mean of the historical deviation divergence set is extracted to construct a historical degradation variation factor. The current health status deviation divergence is then converted into a failure warning probability sequence output by the limit negative exponential decay operator.

[0014] Furthermore, the process of synchronously reading the generator's time-domain active power sequence and time-domain speed sequence under the same sampling clock reference and synchronization timestamp, and obtaining the spatial angular domain equivalent torque sequence with the wind turbine mechanical azimuth angle as the independent variable through continuous integration and spline interpolation resampling, includes:

[0015] The cumulative geometric azimuth is obtained by performing time calculus on the time-domain rotational speed sequence;

[0016] Divide the time-domain active power sequence by the time-domain speed sequence at the same moment to extract the time-domain equivalent torque;

[0017] The conformal spline interpolation operator is invoked to resample and map the time-domain equivalent torque on the time independent variable axis according to the cumulative geometric azimuth angle at equal intervals, thereby generating the spatial angular domain equivalent torque sequence.

[0018] Furthermore, the step of extracting the integral expectation of the equivalent torque sequence in the single-cycle azimuth domain of the wind turbine, with the range of 0 to 2π mechanical angles corresponding to one revolution of the wind turbine as the value range and the azimuth angle corresponding to the Z pulse of the main shaft encoder as the starting point, as the single-cycle average reference torque, and combining it with the third harmonic cosine analytical function to synthesize the third harmonic reference aerodynamic excitation sequence, includes:

[0019] The equivalent torque sequence in the spatial angle domain is calculated by performing calculus in the azimuth domain of a single wind turbine rotation and the average value is obtained to obtain the average reference torque of a single rotation.

[0020] Multiply the absolute geometric azimuth of space by a constant 3, substitute it into the standard cosine analytic function, and take the opposite to generate the reference multiplier;

[0021] The average reference torque per cycle is multiplied and combined with the reference multiplier to synthesize the third harmonic reference aerodynamic excitation sequence.

[0022] Furthermore, the step of performing cross-correlation integration on the third harmonic reference aerodynamic excitation sequence and the spatial angular domain equivalent torque sequence to analyze the torsional hysteresis phase angle of the transmission chain includes:

[0023] A sliding misalignment angle independent variable is set for the azimuth angle domain of a single rotation of the wind turbine;

[0024] The third harmonic reference aerodynamic excitation sequence is multiplied by the spatial angular domain equivalent torque sequence of the sliding misalignment angle independent variable shifted backward, and continuous integration is performed in the corresponding wind turbine single-cycle azimuth domain to generate a continuous cross-correlation function.

[0025] Within the constraint boundary from 0 to 2π / 3, the coordinates of the sliding misalignment angle independent variable corresponding to the maximum value of the continuous cross-correlation function are extracted using the global extremum search operator, and used as the torsional hysteresis phase angle of the transmission chain.

[0026] Furthermore, the ratio of the average integral value of the objective working torque of the current cycle to the torsional hysteresis phase angle of the transmission chain is used as the characteristic value of the dynamic equivalent torsional stiffness in the current period, including:

[0027] Within the azimuth domain of a single wind turbine rotation, the equivalent torque sequence in the spatial azimuth domain is integrally calculated and divided by the full circumference value to obtain the average integral value of the objective working torque of the current cycle.

[0028] Divide the average integral of the objective working torque of the current cycle by the torsional hysteresis phase angle of the transmission chain, and map the extracted deformation hysteresis measure to the current dynamic equivalent torsional stiffness characteristic value.

[0029] Furthermore, the step of solving for the optimal smoothing bandwidth based on the discrete standard deviation of the stiffness feature sample set, constructing the state probability density function, and integrating to generate the system state information entropy includes:

[0030] The mean of the sum of squares of the differences between each current dynamic equivalent torsional stiffness characteristic value and its average value in the stiffness characteristic sample set is calculated, and the discrete standard deviation is obtained by square root operation.

[0031] The optimal smoothing bandwidth is obtained by multiplying the discrete standard deviation by a constant 1.06 and the total number of samples to the power of -1 / 5.

[0032] By combining each stiffness feature sample in the stiffness feature sample set with the optimal smoothing bandwidth, the state probability density function is constructed using a Gaussian kernel distribution;

[0033] The state probability density function is multiplied by a dimensionless corrected constant, and the natural logarithm is taken. This logarithm is then multiplied by the state probability density function itself, and a negative integral is performed over the entire positive half-axis region to generate the system state information entropy.

[0034] Furthermore, the step of retrieving the initial healthy steady-state baseline probability distribution, performing a cross-logarithmic integral transformation on the state probability density function, and deriving the current health state deviation divergence and updating it to the historical deviation divergence set includes:

[0035] The system extracts the baseline probability distribution established under the condition of the equipment being in a non-destructive state during factory trial operation as the initial healthy steady-state reference probability distribution;

[0036] Divide the state probability density function within the current operation window by the initial healthy steady-state baseline probability distribution to obtain the proportional kernel domain;

[0037] After extracting the natural logarithm from the proportional kernel domain, the state probability density function is multiplied by the natural logarithm result. Continuous spatial integration is then performed within the finite ergodic integration domain corresponding to the state probability density function and the initial healthy steady-state baseline probability distribution to calculate the structural performance degradation measure and output the current healthy state deviation divergence, which is then stored in the historical deviation divergence set.

[0038] Furthermore, the step of extracting the squared mean of the historical deviation divergence set to construct a historical degradation variation factor, and then converting the current health status deviation divergence into a failure early warning probability sequence output using a limit negative exponential decay operator, includes:

[0039] The squares of the divergence of health status in each historical divergence set are calculated and the average is summed to construct the historical degradation variation factor;

[0040] Taking the square value of the current health state deviation divergence as the numerator and dividing it by twice the historical degradation variation factor, a negative-valued bounded convergence kernel domain is constructed;

[0041] Substituting the bounded convergence kernel domain as an exponent into an attenuation function with the natural constant as the base, subtracting the output result of the attenuation function from the constant 1, compressing and constraining the continuously increasing divergence measure into the left-closed right-open interval [0,1), and generating and outputting the failure early warning probability sequence.

[0042] The beneficial effects of the present invention are as follows: by resampling and mapping the non-stationary time series to the absolute spatial azimuth coordinates, the spectrum dispersion interference caused by variable rotational speed under natural wind conditions is effectively eliminated; a non-lag base matrix is synthesized by utilizing the unit's own operation expectation and the objective spatial law of the natural flow field, and the hysteresis offset phase angle of mechanical energy transfer is resolved through cross-correlation, and then the actual work load is divided by the phase angle to obtain the dynamic equivalent torsional stiffness characteristic value; this process does not rely on external empirical threshold setting at all, the Gaussian bandwidth and the state probability density function are constructed by extracting the discrete standard deviation, the health state deviation divergence is deduced in combination with the natural thermodynamic entropy theory, and the failure early warning probability sequence is generated by using the historical degradation variation factor, thereby realizing the objective quantification and adaptive early warning of the early fatigue degradation evolution trend inside the unit transmission system without adding additional sensing equipment, and effectively reducing the false alarm rate caused by sudden environmental changes. Brief Description of the Drawings

[0043] Figure 1 is a flowchart of an intelligent monitoring method for wind turbines in a wind farm according to the present invention. Detailed Description of the Embodiments

[0044] The subject matter described herein will now be discussed with reference to example embodiments. It should be understood that the discussion of these embodiments is only to enable those skilled in the art to better understand and implement the subject matter described herein, and the functions and arrangements of the discussed elements can be changed without departing from the protection scope of the content of this specification. Various examples may omit, substitute, or add various processes or components as needed. In addition, features described relative to some examples may also be combined in other examples.

[0045] This embodiment is applied to an electric digital data processing system including a central processing unit. The electric digital data processing system includes an industrial control computer supporting the SCADA system of the wind farm, a cloud computing power cluster, and a local edge computing gateway for the wind turbine. The system realizes bidirectional communication with the main control system of the wind turbine and the main shaft encoder through the industrial Ethernet, the communication protocol adopts ModbusTCP or OPCUA, and the system realizes clock synchronization with all acquisition terminals through the IEEE1588PTP protocol, with a time synchronization error not exceeding 10% of the sampling period.

[0046] like Figure 1 As shown, a method for intelligent monitoring of wind turbines in a wind farm is applied to an electronic digital data processing system containing a central processing unit, comprising:

[0047] Under the same sampling clock reference and synchronization timestamp, the generator time-domain active power sequence and time-domain speed sequence are read synchronously. After continuous integration and spline interpolation resampling, the spatial angular domain equivalent torque sequence with the wind turbine mechanical azimuth angle as the independent variable is obtained.

[0048] Within the azimuth domain of a single wind turbine rotation, with the mechanical angles from 0 to 2π corresponding to one revolution of the wind turbine as the range and the azimuth angle corresponding to the Z pulse of the main shaft encoder as the starting point, the integral expectation of the equivalent torque sequence in the spatial azimuth domain is extracted as the average reference torque for a single revolution, and combined with the third harmonic cosine analytical function to synthesize the third harmonic reference aerodynamic excitation sequence.

[0049] By performing cross-correlation integration on the third harmonic reference aerodynamic excitation sequence and the spatial angular domain equivalent torque sequence, the torsional hysteresis phase angle of the transmission chain is analyzed.

[0050] The ratio of the average integral value of the objective working torque of the current cycle to the torsional hysteresis phase angle of the transmission chain is used as the current dynamic equivalent torsional stiffness characteristic value. The current dynamic equivalent torsional stiffness characteristic value is used to characterize the average torque carrying capacity corresponding to the unit torsional hysteresis phase angle of the transmission chain and is collected into a stiffness characteristic sample set.

[0051] The optimal smoothing bandwidth is calculated based on the discrete standard deviation of the stiffness feature sample set, and the state probability density function is constructed and integrated to generate the system state information entropy.

[0052] Retrieve the initial health steady-state baseline probability distribution, perform a cross-logarithmic integral transformation on the state probability density function, deduce the current health state deviation divergence and update it to the historical deviation divergence set;

[0053] The squared mean of the historical deviation divergence set is extracted to construct a historical degradation variation factor. The current health status deviation divergence is then converted into a failure warning probability sequence output by the limit negative exponential decay operator.

[0054] The system synchronously collects the time-domain active power sequence and time-domain speed sequence from the wind turbine generator side. The sampling period is set to a fixed value in the range of 1ms to 100ms. The sampling period matches the rated speed of the wind turbine. The higher the rated speed, the shorter the sampling period is set. The sampling data format is a single-precision floating-point one-dimensional array with a synchronization timestamp. The active power range is 0 to 1.2 times the rated active power of the unit, and the speed range is 0 to 1.2 times the rated speed of the unit.

[0055] The system executes subsequent calculations only when the wind turbine is in grid-connected steady-state operation. The criteria for grid-connected steady-state operation are: the wind turbine speed is stable between 70% and 105% of its rated speed for at least three rotor rotation cycles; the wind turbine's active power is stable between 20% and 100% of its rated power for at least three rotor rotation cycles; and the wind turbine exhibits no pitch control, no fault alarms, and no unlimited power commands. When the wind turbine is in non-steady-state conditions such as start-up / shutdown, standby, fault ride-through, pitch control, or power-limited operation, the system pauses the calculation process, maintaining the current sample set and historical data unchanged until the operating conditions return to grid-connected steady-state operation.

[0056] The system preprocesses the acquired time-domain active power sequence and time-domain rotational speed sequence. First, it uses the 3σ criterion to identify outliers, classifying sampling points exceeding ±3 standard deviations of the sequence mean as outliers. These outliers are replaced with linear interpolation results from adjacent valid sampling points. Then, a third-order Butterworth low-pass filter is used to filter the sequence, with the filter cutoff frequency set to 10 times the rated rotational frequency of the wind turbine, to remove high-frequency electromagnetic interference and sensor noise. For sequences with no more than 10% of consecutive missing sampling points per revolution, linear interpolation is used to complete the missing data. Sequences with more than 10% of consecutive missing sampling points per revolution are considered invalid data.

[0057] The algorithm uses the single-cycle rotation period of the wind turbine as the smallest execution unit. The determination of the single-cycle rotation period is triggered by the Z-pulse signal from the main shaft encoder. The Z-pulse signal corresponds to the starting zero point of the 0 rad azimuth angle of the wind turbine's single-cycle rotation. Each time a Z-pulse signal is received, the system triggers a complete single-cycle calculation process, completing the calculation of all parameters for the current cycle and updating the corresponding sample set and historical data. In scenarios without a main shaft encoder Z-pulse signal, the starting zero point of the single-cycle rotation is determined by the zero-crossing point of the modulus 2π of the cumulative geometric azimuth angle. The zero-crossing point determination threshold is set to 0.01 rad. When the result of the modulus 2π of the cumulative geometric azimuth angle jumps from greater than 2π - 0.01 rad to less than 0.01 rad, it is determined as the starting zero point of the single-cycle rotation, triggering the single-cycle calculation process.

[0058] The cumulative geometric azimuth is obtained by performing time integration on the time-domain rotational speed sequence. After the system completes the acquisition and preprocessing of the current single-cycle time-domain rotational speed sequence, this step is performed. The calculation formula is as follows:

[0059]

[0060] in, The time corresponding to the starting zero point of a single lap, with the dimension of seconds. The current sampling time is expressed in seconds (s). Let be the time integral variable, with the dimension s. Let the time element be a time infinitesimal with dimensions in seconds. and All are preprocessed time-domain rotational speed sequences, with units of rad / s. The global cumulative geometric azimuth is expressed in rad and represents the cumulative angle of the wind turbine rotation starting from the zero point of a single rotation.

[0061] For discretely sampled time-domain rotational speed sequences, the trapezoidal numerical integration method is used to perform the above continuous integration operation, and the calculation formula is as follows:

[0062]

[0063] in, For the nth sampling time, For the (n-1)th sampling time, The global cumulative geometric azimuth angle corresponding to the nth sampling time. For sampling moments where the rotational speed sequence value is 0, the cumulative geometric azimuth remains unchanged from the previous valid moment. For the end of a single lap, the absolute value of the difference between the modulus 2π of the cumulative geometric azimuth and 2π does not exceed 0.01 rad; otherwise, the single lap data is considered invalid.

[0064] Divide the time-domain active power sequence by the time-domain speed sequence at the same moment to extract the time-domain equivalent torque. After the system completes the preprocessing of the current single-cycle time-domain active power sequence and the calculation of the cumulative geometric azimuth, this step is performed. The calculation formula is:

[0065]

[0066] in, The result is the preprocessed time-domain active power sequence, with dimensions in N·m·rad / s. The time-domain equivalent torque corresponding to the current sampling moment, with dimensions in N·m, characterizes the equivalent driving torque output by the generator side during wind turbine rotation. For wind turbine models with gearboxes, the calculated generator-side time-domain equivalent torque needs to be further multiplied by the gearbox transmission ratio and divided by the generator electromechanical conversion efficiency parameter to align its physical dimensions and map it to the low-speed main shaft side. For sampling moments where the time-domain speed sequence value is less than 0.1% of the rated speed, the corresponding time-domain equivalent torque is directly assigned a value of 0 to avoid division by zero errors. Sampling points where the calculated time-domain equivalent torque exceeds the range of 0 to 1.5 times the rated torque of the unit are determined to be invalid values ​​and replaced with the linear interpolation result of adjacent valid sampling points.

[0067] The conformal spline interpolation operator is invoked to resample and map the time-domain equivalent torque on the time independent variable axis according to the cumulative geometric azimuth angle at equal intervals, generating a spatial angular domain equivalent torque sequence. After the system completes the calculation of the current single-cycle time-domain equivalent torque sequence and the corresponding cumulative geometric azimuth angle sequence, this step is executed. The calculation formula is as follows:

[0068]

[0069] in, The piecewise cubic Hermite interpolation operator is a standard implementation of shape-preserving spline interpolation. It can maintain the monotonicity and concavity of the original data and avoid interpolation spurious oscillations. The angular step size is equidistant and has the dimension of rad. This is a dimensionless step size index, with values ​​ranging from 0 to... -1, This represents the number of equally spaced sampling points per loop. It represents the relative azimuth angle at equal intervals in a single circle, measured in rad, and ranging from 0 to 2π. This is the equivalent torque sequence in the spatial angular domain, with dimensions in N·m, and is the equivalent torque sequence corresponding to the azimuth angles at equal intervals in a single revolution.

[0070] Number of equally spaced sampling points per circle Set to an integer between 360 and 3600. The value of matches the sampling period, satisfying ≥ 0.8 times the original number of sampling points per loop, with equal-interval angular step size. The system uses the cumulative geometric azimuth sequence of the current single cycle as the independent variable and the time-domain equivalent torque sequence as the dependent variable to construct a conformal piecewise cubic Hermite interpolation model. Using the equidistant relative azimuth sequence of the single cycle as the query point, the corresponding spatial angular domain equivalent torque sequence is obtained through interpolation. The boundary conditions for the interpolation calculation are set as periodic boundaries. The derivatives at the start and end positions of a single cycle are obtained by joint central difference calculation using adjacent data points spanning the 0 / 2π start and end points of the wind turbine, ensuring the smoothness and physical continuity of the angular domain signal at the boundary of the rotation cycle. For interpolation results exceeding the reasonable range of the rated torque, the average of adjacent valid interpolation results is used for replacement.

[0071] Within the azimuth domain of a single rotor cycle, the equivalent torque sequence in the spatial angular domain is integrated and its average value is calculated to obtain the average reference torque for that single rotor cycle. After the system generates the equivalent torque sequence in the spatial angular domain for the current single rotor cycle, this step is performed. The average reference torque for that single rotor cycle is the expected integral of the equivalent torque sequence in the spatial angular domain within the current azimuth domain, representing the average aerodynamic load transmitted by the rotor in the current single rotor cycle. For the discrete equivalent torque sequence in the spatial angular domain, the trapezoidal numerical integration method is used to calculate the integral value for that single rotor cycle, which is then divided by the full circumference radius of 2π to obtain the average reference torque for that single rotor cycle.

[0072] Multiply the absolute geometric azimuth angle in space by a constant 3, substitute it into the standard cosine analytic function, and take its opposite to generate the reference multiplier. After the system completes the construction of the single-circle equally spaced relative azimuth angle sequence, this step is performed. The calculation formula is as follows:

[0073]

[0074] in, As a reference multiplier, and a purely proportional dimensionless scalar, the constant 3 corresponds to the third harmonic aerodynamic excitation characteristics brought about by the three-bladed structure of the wind turbine, and matches the periodic load characteristics brought about by wind shear and tower shadow effect during wind turbine rotation. The azimuth angles are equal-spaced relative angles for a single rotation, with dimensions in rad. The phase change of the aerodynamic excitation load during a single rotation of the wind turbine is characterized by the reference multiplier under ideal conditions without torsional hysteresis.

[0075] The single-cycle average reference torque is multiplied and combined with the reference multiplier to synthesize a third-harmonic reference aerodynamic excitation sequence. After the system completes the calculation of the single-cycle average reference torque and the generation of the reference multiplier, this step is performed. The calculation formula is as follows:

[0076]

[0077] in, The angle variable is a single-circle integral, with the dimension rad. Let be an integral element of angle, with dimensions in rad. The average reference torque per revolution is expressed in N·m. The third harmonic reference aerodynamic excitation sequence, with dimensions in N·m, characterizes the third harmonic periodic aerodynamic excitation load sequence borne by the wind turbine under ideal conditions without torsional hysteresis of the transmission chain.

[0078] For discrete sequences, the discrete calculation formula for the average reference torque per lap is:

[0079]

[0080] in, Ns represents the average reference torque per revolution, and Ns represents the number of equally spaced sampling points per revolution. When the average reference torque per revolution is less than 10% of the rated torque, the load for the current revolution is deemed invalid.

[0081] A sliding misalignment angle independent variable is set for the azimuth domain of a single rotor revolution. After the system completes the generation of the third harmonic reference aerodynamic excitation sequence and the equivalent torque sequence in the spatial azimuth domain for the current revolution, this step is performed to calculate the sliding misalignment angle independent variable. The dimension is rad, which represents the relative translation angle between the third harmonic reference aerodynamic excitation sequence and the equivalent torque sequence in the spatial angular domain. It is used to quantify the torque response hysteresis caused by the torsional deformation of the transmission chain. The value ranges from 0 to 2π / 3, and is matched with the period of the third harmonic aerodynamic excitation sequence to avoid confusion of multiple extrema.

[0082] The third-harmonic-frequency reference aerodynamic excitation sequence is multiplied by the spatial angular domain equivalent torque sequence of the backward-translated sliding misalignment angle independent variable, and continuous integration is performed within the corresponding single-rotor azimuth domain to generate a continuous cross-correlation function. After the system completes the setting of the value interval of the sliding misalignment angle independent variable, this step of the calculation is performed. The calculation formula is as follows:

[0083]

[0084] in, This is the index of the current number of rotations in the cycle, dimensionless, and takes the value of a positive integer. The variable is the relative azimuth integral of a single circle, with dimensions in rad and a range of 0 to 2π. Let be an integral element of angle, with dimensions in rad. The spatial angular domain equivalent torque sequence is defined as the independent variable of the backward translational sliding misalignment angle, with dimensions in N·m. Let f(m) be the continuous cross-correlation function corresponding to the m-th cycle, with dimensions N²·m²·rad, representing the strength of the correlation between the two sequences under different translation angles.

[0085] For the translated spatial angular domain equivalent torque sequence, a periodic extension method is used to process the boundary data. hour, ,when hour, For discrete sequences, a uniform grid search method is used, dividing the range of values ​​for the sliding misalignment angle variable into Nsearch search steps, where Nsearch is set to an integer between 100 and 1000. Calculate the cross-correlation function value corresponding to each step size sequentially. The discrete calculation formula is as follows:

[0086]

[0087] in, Let be the sliding misalignment angle corresponding to the i-th search step. The angular step size is equal. The larger the value of the cross-correlation function, the higher the phase matching degree of the two sequences at the corresponding translation angle, and the closer the translation angle is to the actual torsional hysteresis phase angle of the transmission chain.

[0088] Within the constraint boundary range of 0 to 2π / 3, the coordinates of the sliding misalignment angle independent variable corresponding to the maximum value of the continuous cross-correlation function are extracted using the global extremum search operator, and this is used as the torsional hysteresis phase angle of the transmission chain. After the system completes the calculation of the continuous cross-correlation function over the entire interval of the current loop, this step is executed. The calculation formula is as follows:

[0089]

[0090] in, This is a global extremum search operator used to extract the coordinates of the independent variable corresponding to the maximum value of a function. To constrain the boundary, Let be the torsional hysteresis phase angle of the transmission chain corresponding to the m-th turn, with the dimension of rad. It represents the torque response hysteresis angle caused by torsional deformation of the transmission chain under aerodynamic load, and is directly positively correlated with the torsional angle of the transmission chain.

[0091] For discrete cross-correlation function sequences, a local low-frequency detrending process is first performed on continuous cross-correlation functions before extreme value search to filter out macroscopic low-frequency load peaks caused by changes in environmental wind speed. Then, the values ​​corresponding to all search step sizes are traversed, and the sliding misalignment angle corresponding to the maximum value is extracted as the torsional hysteresis phase angle of the transmission chain reflecting microscopic deformation. If multiple identical maximum values ​​exist, the smallest sliding misalignment angle is taken as the result. When the extracted torsional hysteresis phase angle of the transmission chain is less than 1e-5 rad, it is determined to be an invalid value, and the average hysteresis phase angle of the first three valid loops is used to replace it. When the hysteresis phase angle exceeds the constraint boundary of 0 to 2π / 3, the current loop is determined to be invalid and is not included in the subsequent sample set. The calculation in this step is completely consistent with the theory of torsion of circular shafts in classical mechanics of materials. The torsion angle of a circular shaft is directly proportional to the torque and inversely proportional to the torsional stiffness. The torsional deformation of the transmission chain will cause the actual torque response to lag behind the ideal aerodynamic excitation sequence. This lag angle is the torsional hysteresis phase angle, which is equivalent to the torsion angle of the transmission chain and can be used to characterize the torsional stiffness characteristics of the transmission chain.

[0092] Within the azimuth domain of a single wind turbine rotation, the equivalent torque sequence in the spatial angular domain is integrated and divided by the full circumference radian value to obtain the average integral value of the objective working torque for the current cycle. After the system completes the calculation of the equivalent torque sequence in the spatial angular domain and the torsional hysteresis phase angle of the transmission chain for the current cycle, this step is performed. The average integral value of the objective working torque for the current cycle is the same physical quantity as the average reference torque for a single cycle, representing the average driving torque actually output by the wind turbine in the m-th cycle, with dimensions in N·m. The calculation formula is consistent with the average reference torque for a single cycle.

[0093] Divide the average integral value of the objective working torque of the current cycle by the torsional hysteresis phase angle of the transmission chain, and map the extracted deformation hysteresis measure to the current dynamic equivalent torsional stiffness characteristic value. After the system completes the calculation of the average integral value of the objective working torque of the current cycle and the torsional hysteresis phase angle of the transmission chain, this step is performed. The calculation formula is:

[0094]

[0095] in, This represents the integral mean of the objective work torque for the current cycle corresponding to the m-th cycle, with dimensions in N·m. The characteristic value of the dynamic equivalent torsional stiffness corresponding to the m-th cycle is in N·m / rad. It represents the torque carrying capacity of the transmission chain per unit torsional angle. The larger the value, the higher the torsional stiffness of the transmission chain and the smaller the torsional deformation.

[0096] This formula is similar to the classical material mechanics formula for torsional stiffness. Correspondingly, where T is torque, The torsion angle is defined as T, the average torque per revolution is defined as T, and the torsional hysteresis phase angle is defined as T. This enables sensorless real-time calculation of the dynamic torsional stiffness of the transmission chain. When the calculated dynamic equivalent torsional stiffness characteristic value exceeds the range of 30% to 120% of the factory reference stiffness, it is determined to be an invalid value and is not included in the stiffness characteristic sample set. The stiffness characteristic sample set only includes the dynamic equivalent torsional stiffness characteristic value corresponding to the valid loop.

[0097] The stiffness feature sample set is constructed using a sliding time window. The sliding window size is set to 100 to 1000 effective rotations, and the window step size is set to one effective rotation. After each effective rotation, the stiffness feature sample set is updated, removing historical samples outside the window and adding valid samples from the current rotation. The total number of samples in the stiffness feature sample set is denoted as... , dimensionless, and its value is equal to the number of valid samples within the sliding window.

[0098] Calculate the sum of squares and the mean of the differences between each current dynamic equivalent torsional stiffness eigenvalue and its average value in the stiffness feature sample set, and then take the square root to obtain the discrete standard deviation. This step is performed after the system updates the current window's stiffness feature sample set. The calculation formula is:

[0099]

[0100] in, This is an index for stiffness feature samples, dimensionless. Let i be the i-th stiffness feature sample in the stiffness feature sample set, with dimensions in N·m / rad. It is the arithmetic mean of all eigenvalues ​​in the current stiffness feature sample set, with dimensions in N·m / rad. The discrete standard deviation, with dimensions in N·m / rad, characterizes the degree of dispersion of stiffness feature samples.

[0101] Multiply the discrete standard deviation by a constant 1.06 and the total number of samples to the power of -1 / 5 to obtain the optimal smoothing bandwidth. After the system completes the calculation of the discrete standard deviation, this step is performed. The calculation formula is as follows:

[0102]

[0103] in, For optimal smoothing bandwidth, the dimension is N·m / rad, and the constant is 1.06. This is the standard coefficient used to calculate the Silverman optimal bandwidth in Gaussian kernel density estimation. It is used to balance the bias and variance of kernel density estimation and is suitable for data with a unimodal approximate normal distribution. When the stiffness feature sample set exhibits multimodal characteristics, the optimal smoothing bandwidth is multiplied by a correction factor of 0.8 to reduce the bandwidth to adapt to the multimodal distribution and avoid excessive smoothing that could lead to the loss of distribution features.

[0104] Combining the stiffness feature samples within the stiffness feature sample set with the optimal smoothing bandwidth, a state probability density function is constructed using a Gaussian kernel distribution. After the system completes the calculation of the optimal smoothing bandwidth, this step is executed. The calculation formula is as follows:

[0105]

[0106] in, For exponential mapping operators with the natural constant as the base, The stiffness characteristics are traversed through the integration domain variables, with dimensions in N·m / rad. Let be the state probability density function, with dimensions (N·m / rad)⁻¹, representing the probability of the stiffness eigenvalue occurring within different value ranges. The range of the integration domain is set to be from 0.8 times the minimum value to 1.2 times the maximum value of the current stiffness eigenvalue sample set, and the traversal step size is set to 1 / 10 of the optimal smoothing bandwidth to ensure coverage of the distribution range of all valid samples. When the calculated probability density function value is less than 1e-10, it is directly assigned the value of 1e-10 to avoid meaningless results in subsequent logarithmic operations.

[0107] Multiply the state probability density function by a dimensionless corrected constant, take its natural logarithm, and then multiply it by the state probability density function itself. Perform a negative integral over the entire positive half-axis to generate the system state information entropy. After the system has completed the construction of the state probability density function, perform this step. The calculation formula is as follows:

[0108]

[0109] in, This is a dimensionless corrected constant, with a value of 1 N·m / rad, used to eliminate the dimensional influence of the probability density function and ensure the validity of the domain of logarithmic operations. It is an integral infinitesimal element with dimensions N·m / rad. Let be the system state information entropy, a dimensionless constant, representing the dispersion of stiffness feature samples and the uncertainty of the system state. A higher entropy value indicates a greater dispersion of stiffness features and higher uncertainty of the system state. For discrete probability density function sequences, the trapezoidal numerical integration method is used to perform the above continuous integration operation, ensuring that the integration interval matches the range of values ​​traversed within the integration domain.

[0110] The system extracts the baseline probability distribution established under the non-destructive condition of the equipment during factory trial operation as the initial healthy steady-state reference probability distribution. After the system completes the construction of the probability density function of the current window state, it executes this step of the calculation. The initial healthy steady-state reference probability distribution is denoted as... The dimension is (N·m / rad). -1 , which is the probability density distribution of the dynamic equivalent torsional stiffness characteristic value under the ideal healthy state of the equipment without fatigue damage or component failure.

[0111] The initial healthy steady-state baseline probability distribution is constructed using data from the factory trial operation phase after the equipment is installed and commissioned. The trial operation conditions are steady-state grid-connected operation at rated speed and rated power, with a continuous operating time of no less than 72 hours and at least 10,000 valid rotational revolutions. The initial healthy steady-state baseline probability distribution is generated using a construction method completely consistent with the current state probability density function. The ergodic integral domain of the baseline distribution is completely consistent with the ergodic integral domain of the current state probability density function, ensuring a matching support set between the two distributions. The initial healthy steady-state baseline probability distribution is only reconstructed using trial operation data after major component replacements or complete overhauls of the equipment; it is not updated during normal operation.

[0112] The proportional kernel region is obtained by dividing the state probability density function within the current computation window by the initial healthy steady-state baseline probability distribution. After the system completes the retrieval of the initial healthy steady-state baseline probability distribution and aligns the support sets of the two distributions, this step of the computation is performed. The calculation formula is as follows:

[0113]

[0114] in, The proportional kernel is a dimensionless scalar that represents the relative difference between the current state probability density and the baseline healthy state probability density at the same traversal point. For traversal points where the initial healthy steady-state baseline probability distribution value is less than 1e-10, the proportional kernel is directly assigned a value of 1 to avoid division by zero. For traversal points where the current state probability density is non-zero and the ratio increases due to the baseline distribution reaching the 1e-10 lower limit, the upper limit of the proportional kernel is truncated to 100 to represent severe deviation and prevent numerical overflow.

[0115] After extracting the natural logarithm from the proportional kernel domain, the state probability density function is multiplied by the natural logarithm result. A continuous spatial integration operation is then performed within the finite ergodic integration domain corresponding to the common support set of the state probability density function and the initial healthy steady-state baseline probability distribution. This calculates the structural performance degradation measure to output the current healthy state deviation divergence, which is then stored in the historical divergence set. After the system completes the calculation of the proportional kernel domain, this step is performed. The calculation formula is:

[0116]

[0117] in, The current health status deviation divergence is a dimensionless scale used to quantify the degree of difference between the current status distribution and the baseline health status distribution. The higher the value, the greater the deviation between the equipment status and the health baseline, and the more severe the structural performance degradation.

[0118] For discrete sequences, the trapezoidal numerical integration method is used to perform the continuous integration operation described above. The integration interval is consistent with the traversal integration domain of the two distributions. After each calculation of the current health state deviation divergence, it is stored in the historical deviation divergence set in chronological order. The historical deviation divergence set is stored using a lifecycle accumulation method. After the device is powered off and restarted, the historical data is inherited through non-volatile storage media and is not reset. Due to the numerical truncation error in discrete trapezoidal integration, when the calculated current health state deviation divergence has a small negative value, it is forced to zero to ensure that the information entropy is non-negative, and it is stored in the historical deviation divergence set according to the normal process. When the deviation divergence value exceeds 5 times the historical maximum value, it is judged as an outlier and replaced by the average of the first 3 valid deviation divergences.

[0119] The squares of the health status deviation divergences within the historical divergence set are calculated and their mean values ​​are summed to construct a historical degradation variation factor. This step is performed after the system has completed the calculation of the current health status deviation divergence and updated the historical divergence set. The calculation formula is as follows:

[0120]

[0121] in, The total number of valid samples within the historical deviation divergence set, dimensionless. Index for historical deviation divergence samples, dimensionless. Let j be the effective health state deviation divergence within the historical deviation divergence set. The historical degradation variance factor, dimensionless, characterizes the dispersion of health status deviations throughout the equipment's entire lifecycle. The system employs a historical sample inclusion and anti-proliferation mechanism, including deviations only within the range where the probability of a current failure warning is below 0.5. This excludes drastic fluctuations in data from later stages of failure evolution, ensuring that when this variance factor is naturally updated with increasing equipment operating time, it only reflects the normal fluctuations in the equipment's condition due to aging. When the total number of valid samples in the historical deviation variance set is less than 10, the historical degradation variance factor is fixed at 0.1 to avoid calculation distortion caused by insufficient sample size.

[0122] The squared value of the divergence between the current health status and the target health status is used as the numerator, divided by twice the historical degradation variation factor, to construct a bounded convergent kernel region that takes negative values. After the system completes the calculation of the historical degradation variation factor, this step is performed. The calculation formula is as follows:

[0123]

[0124] in, It is a bounded convergent kernel region, a dimensionless scalar, and its value is always less than or equal to 0, ensuring that the result of subsequent exponential operations is within the range of 0 to 1.

[0125] The bounded convergence kernel is substituted as an exponent into a decay function with the natural constant as its base. The output of this decay function is then subtracted from the constant 1 to compress the continuously increasing divergence measure into the left-closed, right-open interval [0,1), generating and outputting a failure warning probability sequence. After the system completes the calculation of the bounded convergence kernel, this step is performed. The calculation formula is as follows:

[0126]

[0127] in, This is a failure warning probability sequence, with values ​​limited to the interval [0,1). It is dimensionless and used to quantify the probability of failure in the equipment's transmission chain. The closer the value is to 1, the higher the risk of failure.

[0128] The system has two levels of early warning thresholds: a first-level threshold of 0.7 and a second-level threshold of 0.95. When the failure warning probability exceeds the first-level threshold for three consecutive calculation windows, the system outputs a first-level warning, prompting maintenance personnel to conduct equipment inspections. When the failure warning probability exceeds the second-level threshold for two consecutive calculation windows, the system outputs a second-level warning, prompting maintenance personnel to shut down the system and conduct component testing. The update frequency of the failure warning probability sequence is consistent with the update frequency of the stiffness feature sample set. The output format is single-precision floating-point numerical values, and the output interface uses an industrial communication protocol compatible with SCADA systems, synchronously outputting to the wind farm monitoring platform.

[0129] This method eliminates spectral dispersion interference caused by natural wind conditions and speed variations by resampling non-stationary time series data onto absolute spatial azimuth coordinates. It synthesizes a hysteresis-free basis matrix using the unit's own operational expectations and the spatial laws of the flow field. By analyzing the hysteresis offset phase angle of mechanical energy transfer through cross-correlation, it obtains the dynamic equivalent torsional stiffness characteristic value. It does not rely on external empirical threshold settings throughout the process. It constructs Gaussian bandwidth and state probability density function through discrete standard deviation, and combines entropy theory to deduce the divergence of health state deviation. It generates a failure early warning probability sequence using historical degradation variation factors. Without adding additional sensing equipment, it realizes the quantification and adaptive early warning of the early fatigue degradation evolution of the machine's transmission system, reducing the false alarm rate caused by sudden environmental changes.

[0130] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.

Claims

1. A method for intelligent monitoring of wind turbines in a wind farm, characterized in that, Applications to electrically digital data processing systems that include a central processing unit, including: Under the same sampling clock reference and synchronization timestamp, the generator time-domain active power sequence and time-domain speed sequence are read synchronously; the cumulative geometric azimuth is obtained by performing time integration on the time-domain speed sequence; the time-domain active power sequence is divided by the time-domain speed sequence at the same moment to extract the time-domain equivalent torque; the conformal spline interpolation operator is called to resample and map the time-domain equivalent torque on the time variable axis according to the cumulative geometric azimuth at equal intervals to generate a spatial angular domain equivalent torque sequence with the wind turbine mechanical azimuth as the independent variable. The conformal spline interpolation operator is a conformal piecewise cubic Hermite interpolation operator. Within the azimuth domain of a single wind turbine rotation, taking the mechanical angles from 0 to 2π corresponding to one revolution of the wind turbine as the range and the azimuth angle corresponding to the Z pulse of the main shaft encoder as the starting point, the integral expectation of the equivalent torque sequence in the spatial azimuth domain is extracted as the average reference torque for a single revolution. The absolute geometric azimuth angle is multiplied by a constant 3, substituted into the standard cosine analytical function, and the opposite is taken to generate a reference multiplier. The average reference torque for a single revolution is multiplied and combined with the reference multiplier to synthesize a third harmonic reference aerodynamic excitation sequence. A sliding misalignment angle independent variable is set for the azimuth domain of a single wind turbine rotation; the third harmonic reference aerodynamic excitation sequence is multiplied by the spatial angular domain equivalent torque sequence that is shifted backward by the sliding misalignment angle independent variable, and continuous integration is performed in the corresponding azimuth domain of the single wind turbine rotation to generate a continuous cross-correlation function; within the constraint boundary from 0 to 2π / 3, the coordinates of the sliding misalignment angle independent variable corresponding to the continuous cross-correlation function reaching its maximum value are extracted using a global extremum search operator, and this coordinate is used as the torsional hysteresis phase angle of the transmission chain; The ratio of the average integral value of the objective working torque of the current cycle to the torsional hysteresis phase angle of the transmission chain is used as the current dynamic equivalent torsional stiffness characteristic value. The current dynamic equivalent torsional stiffness characteristic value is used to characterize the average torque carrying capacity corresponding to the unit torsional hysteresis phase angle of the transmission chain and is collected into a stiffness characteristic sample set. The optimal smoothing bandwidth is calculated based on the discrete standard deviation of the stiffness feature sample set. Combining each stiffness feature sample in the stiffness feature sample set with the optimal smoothing bandwidth, a state probability density function is constructed using a Gaussian kernel distribution. The state probability density function is multiplied by a dimensionless corrected constant, and the natural logarithm is taken. Then, it is multiplied by the state probability density function itself, and a negative integral is performed over the entire positive half-axis domain to generate the system state information entropy. Retrieve the initial health steady-state baseline probability distribution; divide the state probability density function within the current calculation window by the initial health steady-state baseline probability distribution to obtain the proportional kernel domain; extract the natural logarithm from the proportional kernel domain, multiply the state probability density function by the natural logarithm, perform continuous spatial integration within the finite ergodic integration domain corresponding to the common support set of the state probability density function and the initial health steady-state baseline probability distribution, calculate the structural performance degradation measure to output the current health state deviation divergence, and store it in the historical deviation divergence set; The squares of the health state deviation divergences within the historical deviation divergence set are calculated and the average is summed to construct a historical degradation variation factor. The square of the current health state deviation divergence is used as the numerator and divided by twice the historical degradation variation factor to construct a bounded convergence kernel region with negative values. The bounded convergence kernel region is used as an exponent and substituted into a decay function with the natural constant as the base. The output of the decay function is subtracted from the constant 1 to compress and constrain the continuously increasing divergence measure to the left-closed and right-open interval [0,1), generating a failure warning probability sequence and outputting it.

2. The intelligent monitoring method for wind turbines in a wind farm according to claim 1, characterized in that, The step of extracting the integral expectation of the equivalent torque sequence in the single-cycle azimuth domain of the wind turbine, with the range of 0 to 2π mechanical angles corresponding to one revolution of the wind turbine as the value range and the azimuth angle corresponding to the Z pulse of the main shaft encoder as the starting point, as the single-cycle average reference torque, includes: The equivalent torque sequence in the spatial angle domain is calculated by performing calculus on the single-cycle azimuth domain of the wind turbine and the average value is obtained to obtain the single-cycle average reference torque.

3. The intelligent monitoring method for wind turbines in a wind farm according to claim 2, characterized in that, The ratio of the integral average of the current cycle's objective working torque to the torsional hysteresis phase angle of the transmission chain is used as the characteristic value of the current dynamic equivalent torsional stiffness, including: Within the azimuth domain of a single wind turbine rotation, the equivalent torque sequence in the spatial azimuth domain is integrally calculated and divided by the full circumference value to obtain the average integral value of the objective working torque of the current cycle. Divide the average integral of the objective working torque of the current cycle by the torsional hysteresis phase angle of the transmission chain, and map the extracted deformation hysteresis measure to the current dynamic equivalent torsional stiffness characteristic value.

4. The intelligent monitoring method for wind turbines in a wind farm according to claim 3, characterized in that, The step of solving for the optimal smoothing bandwidth based on the discrete standard deviation of the stiffness feature sample set includes: The mean of the sum of squares of the differences between each current dynamic equivalent torsional stiffness characteristic value and its average value in the stiffness characteristic sample set is calculated, and the discrete standard deviation is obtained by square root operation. The optimal smoothing bandwidth is obtained by multiplying the discrete standard deviation by a constant 1.06 and the total number of samples to the power of -1 / 5.

5. The intelligent monitoring method for wind turbines in a wind farm according to claim 4, characterized in that, The baseline probability distribution established by the system under the condition of the equipment being in a non-destructive state during factory trial operation is used as the initial healthy steady-state reference probability distribution.

Images