Fan blade icing state evaluation method based on dynamic weighted mahalanobis distance and gaussian mixture model clustering

Abstract

A method for assessing the icing status of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering includes the following steps: acquiring a wind turbine blade icing fault diagnosis dataset; calculating the dynamic weighted Mahalanobis distance between iced samples and normal samples using normal sample data as a benchmark; mapping the dynamic weighted Mahalanobis distance to a health index of the wind turbine blade's operating status; and using the health index and Gaussian mixture model clustering to quantitatively analyze the icing status of the wind turbine blades, thereby deriving a quantitative classification of the severity of the icing status samples. This method utilizes weighted dynamic Mahalanobis distance to map high-dimensional features to a one-dimensional health index, and then achieves objective classification through unsupervised clustering. It not only provides a new, quantifiable, and interpretable status assessment tool for wind turbine blade icing faults but also offers a new paradigm for predictive health management of complex equipment, possessing significant theoretical and engineering application value.

Description

Technical Field

[0001] This invention relates to the field of wind turbine blade fault assessment technology, specifically to a method for assessing the icing status of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering. Background Technology

[0002] Wind energy, as a highly competitive renewable and clean energy source, demonstrates enormous development potential in wind-rich areas such as high-altitude regions, mountainous areas, and coastal areas. my country boasts vast and widely distributed wind energy reserves and has become the world's largest wind power market. Wind energy is poised to become a crucial energy pillar supporting sustainable development. Wind turbine blades, as one of the core components of wind power generation equipment, are exposed to low temperatures and high humidity during winter, making them more susceptible to icing failures compared to other parts. The impact of icing failures varies depending on the degree of severity: mild icing results in a slight decrease in aerodynamic characteristics and a minor loss in power generation efficiency; moderate icing significantly alters the blade airfoil, reducing lift and output power; severe icing can lead to turbine shutdown, blade structural damage, and even safety accidents. Because wind farms are often located in remote areas, icing failures are difficult to detect and assess in a timely manner.

[0003] In recent years, with the development of technology, data acquisition and monitoring systems for wind turbines have been widely used in the wind power industry and the technology has become relatively mature. Various methods exist for binary classification detection of icing faults on wind turbine blades, utilizing the massive operational status parameters collected in real time by SCADA systems during wind turbine power generation. Among these, data-driven machine learning methods have received widespread attention and application in this field due to their excellent feature extraction and representation capabilities, as well as their end-to-end detection convenience. However, most existing research methods only achieve binary classification detection of "iced" and "not iced," failing to quantify and grade the severity of icing. In data-driven models, considering the severity of icing can better characterize the health of wind turbine blades, providing strong theoretical support for developing differentiated de-icing strategies and implementing preventative maintenance plans.

[0004] Therefore, developing an efficient quantitative assessment technology for the degree of icing faults in wind turbine blades has significant engineering application value and scientific significance for ensuring the safe and stable operation of wind power systems, improving power generation efficiency, and achieving targeted maintenance. Summary of the Invention

[0005] This invention proposes a method for assessing the icing status of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering. This method uses weighted dynamic Mahalanobis distance to map high-dimensional features into a one-dimensional health index, and then achieves objective classification through unsupervised clustering. This not only provides a new quantifiable and interpretable state assessment tool for wind turbine blade icing failures, but also provides a new paradigm for predictive health management of complex equipment, which has important theoretical significance and engineering application value.

[0006] The technical solution adopted in this invention is as follows:

[0007] A method for assessing the icing status of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering includes the following steps: Step 1: Obtain the wind turbine blade icing fault diagnosis dataset; Step 2: Using normal sample data as a benchmark, calculate the dynamic weighted Mahalanobis distance between the iced samples and the normal samples; Step 3: Map the dynamically weighted Mahalanobis distance to a health index of the wind turbine blade operating status; Step 4: Based on the health index and Gaussian mixture model clustering, the health index (HI) of the icing state of wind turbine blades is quantitatively analyzed, and then the severity level of the icing state samples is quantitatively classified.

[0008] In step 1, the wind turbine blade icing fault diagnosis dataset is obtained; specifically as follows: The wind turbine blade icing fault diagnosis dataset refers to raw SCADA monitoring data that has been preprocessed (cleaned, labeled, and balanced) and then input into a machine learning diagnostic model for training and testing. This invention uses Gradient Boosting Decision Tree (GBDT). The diagnostic model outputs a complete dataset consisting of normal samples and icy samples, along with their corresponding feature vectors. For example... Figure 2 The diagram shows a flowchart for obtaining a dataset for wind turbine blade icing fault diagnosis.

[0009] The original dataset for wind turbine blade icing was collected from SCADA system monitoring records of a real wind farm from November to December of a certain year, covering the time periods corresponding to normal operation and icing status of the units. The monitored variables comprehensively cover the operating status of all key subsystems of the wind turbine, including numerical characteristic parameters closely related to blade icing issues, such as temperature difference, yaw position, power, and average blade angle. The wind speed-to-power ratio, average pitch motor temperature, ambient wind speed squared, engine speed, nacelle temperature, average blade DC charging power supply temperature, and average blade charger DC current are among the twelve characteristics listed in Table 1. These variables reflect the actual operating conditions of the unit from multiple dimensions.

[0010]

[0011] The samples from unlabeled periods in the original wind turbine blade icing dataset are not discussed. This means that during the transition from normal operation to icing, there may be a sample interval where the operating condition is difficult to define clearly. This data is usually introduced by system shutdowns or human intervention, exhibiting obvious atypical operating characteristics; it is removed to improve the reliability of model training. On the other hand, box plotting is used to identify and remove noisy data. Furthermore, in actual operation, wind turbines are mostly in a non-icing normal operating state, with blade icing only occurring occasionally. In this case, the ratio of normal operation samples to icing samples is approximately 15:1, a typical case of class imbalance. After data balancing, the ratio of normal to icing samples is adjusted to approximately 1:1, thereby optimizing the distribution structure of the model training data. The preprocessed original wind turbine blade icing dataset is then input into the machine learning model for fault diagnosis, resulting in the wind turbine blade icing fault diagnosis dataset. The model outputs a diagnostic label (normal / iced) for each sample. The machine learning model used in this case is a gradient boosting decision tree model. The final structure of the wind turbine blade icing fault diagnosis dataset is shown in Table 2:

[0012] Note: For ease of display, only some features are listed in Table 2; all 12 features are actually included.

[0013] In step 2, a normal sample refers to SCADA operation data collected when the surface of the wind turbine blades is free of ice. At this time, the aerodynamic shape of the blades is intact and the unit's operating parameters are within the healthy baseline range.

[0014] Icing samples refer to data collected when the surface of wind turbine blades is covered with ice. In this case, the aerodynamic performance of the blades is reduced and the operating parameters deviate from the normal mode. The raw dataset of wind turbine blade icing will have a label column (such as label or icing_status): 0 indicates normal samples and 1 indicates icing samples.

[0015] Step 2 includes: 2.1: Calculation of weighted Mahalanobis distance: Information entropy is used to evaluate the information content of each feature indicator, thereby characterizing its effectiveness in practical applications. Information entropy can be calculated based on the degree of variation of each feature dimension, assigning it its entropy value and corresponding weight. n Includes p For samples with 3D features, the steps for calculating information entropy and weights are as follows: 1): For the first i The first sample j Features Data standardization is performed, and the data is compressed into numbers within the range [0, 1] using the Z-score method: (1); (2); (3); In the above formula: Representing the i The sample at the th j Features The standard value above; Indicates the first j The mean of a feature across all samples. Indicates the first j The standard deviation of each feature. This represents the total number of samples.

[0016] 2): Calculate the first j The first feature i The proportion of each sample, i.e., the probability matrix for calculating relative entropy: (4); In equation (4): Indicates the first j The first feature i The proportion of each sample.

[0017] 3): Calculate the information entropy and information utility value of each feature, and then normalize to obtain the entropy weight of each feature: For the j For each feature, the formula for calculating its information entropy is: (5); In equation (5): Indicates the first j Information entropy of each feature; Represents the probability logarithm with base natural logarithm, when hour, .

[0018] The calculation formula is: (6); In formula (6): Indicates the first j Features .

[0019] After normalization, the entropy weights for each feature are obtained, as shown in the following formula: (7); In equation (7): Indicates the normalized i-th j Features ; This represents the total number of features.

[0020] 4) Furthermore, the weight of a feature in the model should match its sensitivity to abnormal states. That is, features containing more variation information are more beneficial for state identification and prediction. Therefore, Fisher discriminant analysis is used to determine the feature sensitivity weights. for: (8); In equation (8): Indicates the first The sensitivity of each feature and The first and second samples were respectively normal samples and iced samples. The mean of each feature; and These represent the total number of normal samples and the total number of icing samples, respectively. and The first and second samples were respectively normal samples and iced samples. Observations of each feature; The meaning of the first A normal sample, Indicates the first j One indicator, Indicates the first One ice-covered sample.

[0021] Fisher's criterion function value As a metric, the decision weights for calculating features are: (9); In equation (9): Represents the first term calculated using Fisher's criterion. The decision weights of each feature; Indicates the first The sensitivity of each feature; Indicates the first One feature; This represents the total number of features.

[0022] 5): Combining information entropy and Fisher's criterion, the comprehensive weight of each attribute is calculated as follows: (10); In formula (10): Indicates the first The combined weight of each feature; The first term calculated using information entropy is represented by... The information weights of each feature.

[0023] The weighted Mahalanobis distance for all samples in the normal class is defined as: (11); In equation (11): Indicates the first Weighted Mahalanobis distance of normal samples This indicates the first in the normal class. Standardized data for each sample; This represents the standard deviation of a normal sample. This represents the matrix transpose operation.

[0024] ) is the feature weight matrix of the normal class samples; They represent the first The comprehensive weight of the first feature, the second feature The combined weight of each feature, No. The combined weight of each feature.

[0025] (12); In equation (12): This indicates the total number of ice-covered samples.

[0026] (13); The weighted Mahalanobis distance for all samples in the corresponding icing class is: (14); In equation (14): Indicates the first Weighted Mahalanobis distance for each ice-covered sample. Indicates the first in the ice-covered category Standardized data for each sample.

[0027] 2.2: The weighted dynamic Mahalanobis distance is calculated as follows: Assumption for A normal sample, in, express ; These represent the first normal sample and the second normal sample, respectively. .

[0028] Z These are standardized, normal sample data.

[0029] This represents a vector augmented matrix.

[0030] The augmented matrix of a normal sample can be expressed as follows: , in, This represents the total augmented matrix at different lag times for normal samples; Let represent the augmented matrix at the first lag time of the normal sample, and represent the augmented matrix at the second lag time of the normal sample, respectively. Normal sample number An augmented matrix with a time lag; q This indicates the maximum lag order.

[0031] Therefore, the weighted dynamic Mahalanobis distance (WDMD) for normal samples is: (15); In equation (15): Indicates the first Weighted dynamic Mahalanobis distance of a normal sample This represents the standard deviation of the normal sample.

[0032] (16); In equation (16): Indicates the first in the normal class Standardized data for each sample. This represents the total number of normal samples.

[0033] (17); In equation (17): This indicates the matrix transpose operation.

[0034] Let the augmented matrix of the icing sample be... , in: This represents the total augmented matrix at different lag times for the icing sample; Let represent the augmented matrix at the first lag time of the icing sample, and represent the augmented matrix at the second lag time of the icing sample, respectively. , frozen sample An augmented matrix with a lag time.

[0035] The weighted dynamic Mahalanobis distance (WDMD) for the icing samples is calculated as follows: (18); In equation (18): Indicates the first Weighted dynamic Mahalanobis distance for normal samples.

[0036] (19); In equation (19): This represents the standard deviation of the normal sample. This represents the total number of normal samples. Indicates the first in the normal class Standardized data for each sample.

[0037] (20); In equation (20): Indicates the first in the normal class Standardized data for each sample.

[0038] In step 3, the weighted dynamic Mahalanobis distance is mapped to the health status index of the wind turbine blades, and the calculation is as follows: 3.1: Design of Icing Index for Wind Turbine Blades Based on Weighted Dynamic Mahalanobis Distance (WDMD) As shown in the following formula: (twenty one); In equation (21): Represented by natural constant e As the base, with and The product is an exponentiation operation.

[0039] This indicates a regulatory index; Representing all samples The weighted dynamic Mahalanobis distance at time; because ,so ∈(0,1); (twenty two); In equation (22): This represents the confidence level under the condition that the wind turbine blades are healthy, normal, and not iced. The average of the normal samples WDMD ; (twenty three); In equation (23): The average of normal samples WDMD .

[0040] A clear mapping relationship is established between the constructed HI and WDMD, which can effectively reflect the evolution trend of the system state.

[0041] 3.2: The WDMD of normal samples generally does not follow a normal distribution. In order to determine... To obtain the WDMD value, the WDMD of normal samples must first be transformed into a normal distribution using appropriate transformations. The process of processing WDMD values ​​using the Box-Cox transformation is shown in the formula. The transformation formula is: (twenty four); In equation (24): The WDMD values ​​of normal samples after Box-Cox transformation. This represents a power transform parameter (also called a transform coefficient). It is estimated by the maximum likelihood method; First, construct the likelihood function. : (25); In equation (25): yes Maximum likelihood estimation; Indicates a given When dealing with sample data, the likelihood value of the transformed sequence following a normal distribution. This represents the WDMD value of a normal sample before the Box-Cox transformation.

[0042] For all of equation (25) ,have: (26); In equation (26): Indicates the first one sample W MD Take (λ) 1) Power; This represents the total number of normal samples.

[0043] For each In Equation (26) Equation (27) can be expressed as: (27); In equation (27): Normal samples after Box-Cox transformation The average value; (28); In equation (28): This represents the total number of normal samples. Indicates the first in the normal class Standardized data for each sample. Indicates the index of a normal sample.

[0044] The following equation can be derived: (29); Before conversion It does not conform to a normal distribution. After Box-Cox transformation... The distribution characteristics satisfy the requirements of a normal distribution. Using the Box-Cox transform and... Criteria for determining thresholds As shown in equation (30), the confidence level of the wind turbine blades under healthy, normal, and non-icing conditions can then be determined. value.

[0045] (30); In equation (30): Normal samples after Box-Cox transformation The average value; For normal samples after Box-Cox transformation Standard deviation.

[0046] In step 4, a health status index threshold is calculated based on Gaussian mixture model (GMM) clustering, and then the health status index of the wind turbine blade icing state is quantitatively analyzed to obtain a quantitative classification of the severity of the icing state samples, specifically including: S4.1: Input data: Obtain the health status index and label reflecting the health level of the equipment from the sample: 0 indicates normal, 1 indicates icing; S4.2: Determine the optimal number of Gaussian clusters: Use the Bayesian Information Criterion (BIC) as the basis for model selection, as shown in the following equation (31): (31); In equation (31): BIC(k) This represents the number of clusters determined based on the Bayesian information criterion. L It is the maximum value of the likelihood function corresponding to the model. K It is the number of clusters. n It refers to the number of samples.

[0047] As shown in equation (31), the higher the goodness of fit or the smaller the penalty value, the smaller the BIC value. Therefore, the model with the smallest BIC value is selected as the optimal model. The optimal number of clusters is selected as the number of clusters with the smallest BIC value: (32); In equation (32): , Both represent the optimal number of clusters.

[0048] S4.3: GMM Training: Using the Optimal Number of Clusters Fit and train the GMM, then predict the cluster label and probability of each sample, and obtain the cluster label of each sample and the probability of each sample belonging to each Gaussian component: Initialize mean Covariance Weight In the iteration, each observation is calculated. Belonging to the k The posterior probability (also called liability) of a probability distribution function is: (33); In equation (33): Indicates the first i Sample Belonging to the k The posterior probability (also known as liability) of each Gaussian component. For the first k The probability density function of a Gaussian distribution with Gaussian components; For the first k The mean of the distribution of the Gaussian components; For the first k Covariance of Gaussian components; For the first k The weights of each Gaussian component; This represents the total number of clusters.

[0049] (34); In equation (34): Let represent the Gaussian kernel function, and represent the sample. With the k The center of each Gaussian component Similarity measure between them; Indicates the first i Observations of a sample; Indicates the first k The mean of the Gaussian components.

[0050] Then, the following parameter values ​​are calculated: (35); (36); (37); S4.4: Obtain the mean of the GMM Standard deviation and weight : (38); In equation (38): They represent the first 1The mean of the nth Gaussian components, the th 2 The mean of the Gaussian components, , No. k The mean of the Gaussian components.

[0051] (39); In equation (39): They represent the first 1 The covariance of the Gaussian components, the first 2 The covariance of the Gaussian components, , No. k The covariance of the Gaussian components.

[0052] (40); In equation (40): They represent the first 1 The weight of the Gaussian component, the first 2 The weights of the Gaussian components, , No. k The weights of each Gaussian component.

[0053]

[0054] In equation (41): the Gaussian mixture model is used to obtain K One Gaussian component, parameters calculated by mean. Sort in ascending order, the corresponding standard deviation is The weight is ,and =1.

[0055]

[0056] In equation (42): for two adjacent Gaussian components and +1, we need to find a threshold. t This makes the weighted probability densities equal, where... Let be the probability density function of a Gaussian distribution (normal distribution). .

[0057] =

[0058] In equation (43): Represents the Gaussian kernel function (exponential part), and represents the threshold. With the k The center of each Gaussian component Similarity measure between Indicates the first k The mean (center position) of the Gaussian components. For the first k The covariance of the Gaussian components.

[0059]

[0060] In equation (44): Taking the natural logarithm of both sides of equation (43), we can rearrange it into a quadratic equation. a= (45); In equation (45): a Indicates the first k The Gaussian component and the first k+ The difference in variance of one Gaussian component. For the first k The standard deviation of each Gaussian component.

[0061] b= (46); In equation (46): b Indicates the first k The Gaussian component and the first k+ The difference in the combined ratio of the mean to the variance of a single Gaussian component. Indicates the first k+ The mean of a single Gaussian component. Indicates the first k+ The standard deviation of one Gaussian component.

[0062] c= (47); In equation (47): c This represents the effect of the combined weights and standard deviations of the two components on the location of the intersection point. Indicates the first k The weights of each Gaussian component. Indicates the first k+ The weight of a Gaussian component.

[0063] Discriminant: (48); In equation (48): The discriminant in the quadratic equation for threshold calculation in expression (44).

[0064] if ≥0: (49); (50); In equations (49) and (50), This represents one of the solutions to the quadratic equation in expression (44). The expression (44) represents another solution to the quadratic equation.

[0065] (51); In equation (51), Indicates in and Choose a solution from among them.

[0066] S4.6: Based on the calculated threshold, the icing health index is used as the quantitative standard to realize the quantitative classification of the icing health index of wind turbine blades, thereby dividing the samples into different severity levels. According to the threshold, they are classified from low to high as ['severe icing', 'moderate icing', 'mild icing', 'minor icing'].

[0067] (52); S4.7: Count the number of samples for each severity level: (53); In formula (53): This represents the sample size for each severity level. This indicates the number of samples for each severity level. Indicates a statistical quantity.

[0068] This invention provides a method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering. The technical advantages are as follows: 1) Step 2 of this invention proposes a method for constructing a weighted dynamic Mahalanobis distance (WDMD), the core advantages of which include: ① Multi-dimensional Feature Adaptive Weighting: Unlike traditional Mahalanobis distance which applies equal weight to each feature, this step integrates information entropy weights and Fisher's discriminative weights. Information entropy weights are automatically assigned based on the inherent variability (information content) of each feature, avoiding the subjectivity of manually setting weights. Fisher's weights enhance the contribution of icing-sensitive features based on the inter-class separability between normal and iced samples. The combined weights reflect both the intrinsic information content of the features and highlight key indicators for distinguishing icing states, significantly improving the discriminative power of the distance metric.

[0069] ② Introducing dynamic timing characteristics: Traditional Mahalanobis distance only considers static samples at a single time moment, failing to capture the cumulative effect of the icing process over time. This step involves constructing an augmented matrix (introducing the maximum lag order). q This method incorporates the feature vectors from the current moment and several previous moments into the calculation, forming a weighted dynamic Mahalanobis distance. It can reflect the continuous deterioration process of icing conditions, rather than detecting anomalies at isolated moments.

[0070] It exhibits higher sensitivity and robustness to slowly changing faults (such as blade icing).

[0071] ③Standardized benchmarks and comparability: Using the mean and covariance matrix of normal samples as benchmarks, the WDMD of all icing samples is calculated relative to the health benchmark, making the distance values ​​between different operating conditions and different wind turbines horizontally comparable, providing a unified quantitative scale for the subsequent construction of the health index.

[0072] ④ Reduce noise and outlier interference: The weighted Mahalanobis distance itself has the advantage of considering feature correlation. Combined with dynamic expansion, it can effectively suppress the drastic impact of instantaneous noise on the distance value, making the WDMD curve smoother and more trend-oriented.

[0073] 2) Step 3 of this invention maps WDMD to the Health Index (HI), which has the following advantages: ① Nonlinear bounded mapping: using the Sigmoid function HI WDMD is compressed from [0,+∞) to the (0,1] interval. When WDMD is very small (healthy state), HI is close to 1; when WDMD is very large (severe icing), HI approaches 0. The mapping result intuitively reflects the "health level", which is convenient for engineers to understand and make decisions.

[0074] ② Dynamic adjustment parameters α Physical meaning of parameters: α Confidence level of normal samples β The relationship with the average WDMD of normal samples was determined. This ensures that under normal conditions... HI The HI value falls in the high value range (e.g., above 0.95), while the HI value of the icy samples decreases significantly.

[0075] The steepness of the mapping curve is adaptively adjusted to accommodate the sensitivity differences of different wind farms or turbines.

[0076] ③ Box-Cox Transform Ensures the Reliability of Statistical Thresholds: WDMD in normal samples often does not follow a normal distribution, and directly calculating quantiles introduces bias. This step uses the Box-Cox transformation to convert WDMD into an approximately normal distribution, and then uses the 2σ criterion to determine the statistical boundaries of the health status. This improves the scientific rigor and accuracy of setting the health baseline threshold. The transformed data satisfies the assumptions of parametric statistical methods, making it easy to extend to other statistical control chart methods.

[0077] ④ The health index has both clear physical and statistical significance: HIIt reflects both the magnitude of the Mahalanobis distance of the sample relative to the healthy baseline and provides the probabilistic significance of exceeding the normal range through statistical transformation (e.g., 99.7% confidence). This dual interpretability makes the assessment results more credible.

[0078] 3) Step 4 of this invention uses a Gaussian mixture model (GMM) to perform unsupervised clustering of the health index, achieving an objective classification of the severity of icing. The main advantages are as follows: ① Unsupervised automatic classification avoids subjective threshold division: Traditional methods often rely on expert experience or fixed percentages to classify severity levels, lacking data-driven objectivity. This step utilizes Gaussian Mixture Models (GMMs) to perform soft clustering of the HI distribution, automatically identifying the inherent multimodal structure of the data and classifying it into four levels: "mild, moderate, moderate, and severe." The clustering results are entirely determined by the data's own distribution characteristics, making them reproducible and unbiased.

[0079] ② Bayesian Information Criterion (BIC) automatically determines the optimal number of clusters: by calculating different numbers of clusters... k The optimal BIC value is selected by choosing the model with the smallest BIC. This avoids the arbitrariness of manually specifying the number of clusters and also prevents overfitting or underfitting.

[0080] ③ Thresholds are calculated based on the intersection points of Gaussian components, which has probability balance: Instead of using a simple mean midpoint, this method finds the intersection point where the weighted probability density functions of two adjacent Gaussian components are equal. This intersection point corresponds to the decision boundary with the minimum misclassification probability for the two classes of samples, which is mathematically optimal. When the two components have unequal variances or different weights, the threshold automatically biases towards the component with smaller variance or larger weight, which better reflects the actual data distribution.

[0081] ④ Output probabilistic membership degrees to support uncertainty representation: GMM not only provides the clustering label for each sample, but also outputs its posterior probability (responsibility value) of belonging to each Gaussian component. For samples with ambiguous boundaries, it can be determined to what extent they belong to a certain level, providing richer information for subsequent risk decisions.

[0082] ⑤ Count the number of samples for each severity level: The final output is the sample count for each severity level, which can be directly used for wind farm icing event statistics, operation and maintenance plan formulation, and de-icing strategy optimization.

[0083] 4) This invention addresses the limitation of most models and methods that only achieve binary classification detection of "iced" and "uniced," failing to quantify and grade the severity of icing. To achieve a quantitative assessment of icing severity, this invention inputs a wind turbine blade icing fault diagnosis dataset, constructs a weighted dynamic Mahalanobis space using normal samples from the dataset, maps high-dimensional features to a one-dimensional health index using the weighted dynamic Mahalanobis distance, calculates the health status index using WDMD, and then uses the GMM model for unsupervised clustering to achieve icing severity grading (mild). Mild moderate (Severe), providing a quantifiable assessment tool for operation and maintenance decisions, with a novel approach and engineering practicality. This breaks through the current research's predominantly binary diagnostic approach, paving the way for differentiated solutions. Precise maintenance decisions provide a direct basis for decision-making, directly serving the preventative maintenance and intelligent operation and maintenance of wind farms. Through accurate icing assessment, optimal de-icing strategies can be formulated, such as activating them only when necessary, effectively reducing power generation losses. Studies show that precise monitoring can help reduce losses by 8%-12%, lower the energy consumption and equipment wear of ineffective de-icing, ensure wind turbine safety, and improve the economic benefits of the wind farm throughout its entire lifecycle. Attached Figure Description

[0084] The present invention will be further described below with reference to the accompanying drawings and examples; Figure 1 The flowchart of the method of the present invention Figure 2 A flowchart for obtaining a dataset for wind turbine blade icing fault diagnosis.

[0085] Figure 3(a) shows the MD plot of normal samples in the test set; Figure 3(b) shows the MD plot of the icing samples in the test set.

[0086] Figure 4(a) shows the WMD plot of normal samples in the test set; Figure 4(b) shows the WMD plot of the icing samples in the test set.

[0087] Figure 5(a) shows the WDMD plot of normal samples in the test set with lag order q=1; Figure 5(b) shows the WDMD plot of the icing samples in the test set with hysteresis order q=1.

[0088] Figure 6(a) shows the lag order. q =2 WDMD plots of normal samples in the test set; Figure 6(b) shows the lag order. q =2 WDMD plots of icing samples in the test set.

[0089] Figure 7(a) shows the lag order. q=WDMD plots of normal samples in the test set; Figure 7(b) shows the lag order. q =3 WDMD plots of icing samples in the test set.

[0090] Figure 8(a) shows the lag order. q =4 WDMD plots of normal samples in the test set; Figure 8(b) shows the lag order. q =4 WDMD plots of icing samples in the test set.

[0091] Figure 9(a) shows the lag order. q =5 WDMD plots of normal samples in the test set; Figure 9(b) shows the lag order. q =5 WDMD plots of icy samples in the test set.

[0092] Figure 10(a) shows the lag order. q =WDMD plots of normal samples in the test set; Figure 10(b) shows the lag order. q =6 WDMD plots of icy samples in the test set.

[0093] Figure 11 This is a normal distribution plot of WDMD before conversion.

[0094] Figure 12 This is the normal distribution plot of WDMD after conversion.

[0095] Figure 13 The optimal number of clusters for GMM clustering determined by the BIC criterion.

[0096] Figure 14 A visualization analysis chart of the HI quantification results for the icing state of wind turbine blades in the test set. Detailed Implementation

[0097] A method for assessing the icing status of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering. First, normal samples are extracted from the wind turbine blade icing fault diagnosis dataset. A weighted dynamic Mahalanobis space is constructed using these normal samples. The weighted dynamic Mahalanobis distance between iced samples and normal samples is calculated. A health index model is then constructed using the weighted dynamic Mahalanobis distance. Finally, unsupervised clustering using the Gaussian mixture model is applied to classify the icing severity: mild. Mild moderate This approach provides a quantifiable assessment tool for operational and maintenance decisions, offering a novel approach with practical engineering applications.

[0098] Example: Compared to traditional MD WMD The performance of the WDMD test set state evaluation method is shown in Table 3, and Figures 3(a), 3(b) to 10(a), and 10(b). The comparison results in the figures show that when the lag order is set to... q When = 2, the WDMD method performs best in distinguishing between normal and icing states of wind turbine blades and has high computational efficiency. In contrast, the WDMD values ​​of normal and icing samples obtained by traditional MD or WMD methods have a large degree of overlap, making it difficult to effectively separate the two states. To further quantify the distinguishing ability of each method, an overlap rate calculation formula can be introduced for evaluation, see formula (54).

[0099] (54); In the formula: For in the range The number of normal samples in the test set. For in the range The number of icing samples in the test set. For in the range The number of normal samples in the test set. For in the range The number of icing samples in the test set.

[0100] It can be seen that the WDMD calculated according to the method of the present invention ( q≥2 (At the time), the overlap rate is significantly lower than that of traditional MD and WMD, proving that the method proposed in this invention can clearly distinguish between the normal and icing states of wind turbine blades, with a differentiation improvement rate of 23.65% compared to traditional MD and 22.96% compared to WMD.

[0101]

[0102] Before conversion It does not conform to a normal distribution, and after Box-Cox transformation The distribution characteristics satisfy the requirements of a normal distribution. For example... Figure 11 , Figure 12 As shown, and >0, according to experiments =-1.373.

[0103] The optimal number of Gaussian clusters is selected using the BIC criterion. Figure 13 As shown. By Figure 13 It can be seen that the optimal number of clusters k =5. Based on the threshold calculated by GMM clustering, the HI quantification level of wind turbine blade icing status was achieved, as shown in Table 4, thereby dividing the samples into different severity levels.

[0104]

[0105] The sample size for each severity level is shown in Table 5:

[0106] Figure 14 The visualization of icing condition assessment results includes the GMM probability density function and various icing degrees. HI Value distribution, sample classification scatter plot, severity statistics pie chart, such as Figure 14 The research findings not only provide a theoretical basis and technical support for preventive maintenance and intelligent operation and maintenance of wind turbine icing, but also provide a new paradigm that integrates mechanisms and data for multi-state health assessment of complex equipment.

Claims

1. A method for assessing the icing status of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering, characterized in that... Includes the following steps: Step 1: Obtain the wind turbine blade icing fault diagnosis dataset; Step 2: Using normal sample data as a benchmark, calculate the dynamic weighted Mahalanobis distance between the iced samples and the normal samples; Step 3: Map the dynamically weighted Mahalanobis distance to a health index of the wind turbine blade operating status; Step 4: Based on the health index and Gaussian mixture model clustering, a quantitative analysis of the icing status of wind turbine blades is achieved, thereby obtaining a quantitative classification of the severity of the icing status samples.

2. The method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering as described in claim 1, characterized in that: In step 1, the raw dataset of wind turbine blade icing is collected from the SCADA system monitoring records of the wind farm, covering the time periods corresponding to normal operation and icing status of the unit; the monitored variables comprehensively cover the operating status of each key subsystem of the wind turbine unit, including numerical characteristic parameters closely related to the blade icing problem: temperature difference, yaw position, power, average angle of wind turbine blades, etc. Wind speed-to-power ratio, average temperature of pitch motor, square of ambient wind speed, engine speed, nacelle temperature, average temperature of blade DC charging power supply, and average DC current of blade charger. The preprocessed raw dataset of wind turbine blade icing is input into a machine learning model for fault diagnosis, which yields a wind turbine blade icing fault diagnosis dataset. The model outputs the diagnostic label for each sample: normal or icing.

3. The method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering as described in claim 2, characterized in that: In step 2, a normal sample refers to SCADA operation data collected when the surface of the wind turbine blades is free of ice. At this time, the aerodynamic shape of the blades is intact and the unit's operating parameters are within the healthy baseline range. Icing samples refer to data collected when the surface of the wind turbine blades is covered with ice, at which point the aerodynamic performance of the blades decreases and the operating parameters deviate from the normal mode. The original dataset of wind turbine blade icing will have a label column: 0 represents normal samples and 1 represents icing samples.

4. The method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering as described in claim 3, characterized in that: Step 2 includes: calculating the weighted Mahalanobis distance. Information entropy is used to evaluate the information content of each feature indicator to characterize its effectiveness in practical applications; information entropy can be calculated based on the degree of variation of each feature dimension, along with its corresponding weight; for n Includes p For samples with 3D features, the steps for calculating information entropy and weights are as follows: 1): For the first i The first sample j Features Data standardization is performed, and the data is compressed into numbers within the range [0, 1] using the Z-score method: （1）； （2）； （3）； In the above formula: Representing the i The sample at the th j Features The standard value above; Indicates the first j The mean of each feature across all samples; Indicates the first j The standard deviation of each feature; Indicates the total number of samples; 2): Calculate the first j The first feature i The proportion of each sample, i.e., the probability matrix for calculating relative entropy: （4）； In equation (4): Indicates the first j The first feature i The proportion of each sample; 3): Calculate the information entropy and information utility value of each feature, and then normalize to obtain the entropy weight of each feature: For the j For each feature, the formula for calculating its information entropy is: （5）； In equation (5): Indicates the first j Information entropy of each feature; Represents the probability logarithm with base natural logarithm, when hour, ; The calculation formula is: （6）； In formula (6): Indicates the first j Features ; After normalization, the entropy weights for each feature are obtained, as shown in the following formula: （7）； In equation (7): Indicates the normalized i-th j Features ; Indicates the total number of features; 4) Furthermore, the weight of a feature in the model is matched to its sensitivity to abnormal states; that is, the more variation information a feature contains, the better it is for state identification and prediction; feature sensitivity weights are determined based on Fisher discriminant analysis. for: （8）； In equation (8): Indicates the first The sensitivity of each feature and The first and second samples were respectively normal samples and iced samples. The mean of each feature; and These represent the total number of normal samples and the total number of icing samples, respectively. and The first and second samples were respectively normal samples and iced samples. Observations of each feature; The meaning of the first A normal sample, Indicates the first j One indicator, Indicates the first One ice-covered sample; Fisher's criterion function value As a metric, the decision weights for calculating features are: （9）； In equation (9): Represents the first term calculated using Fisher's criterion. The decision weights of each feature; Indicates the first The sensitivity of each feature; Indicates the first One feature; Indicates the total number of features; 5): Combining information entropy and Fisher's criterion, the comprehensive weight of each attribute is calculated as follows: （10）； In formula (10): Indicates the first The combined weight of each feature; The first term calculated using information entropy is represented by... Information weights of each feature; The weighted Mahalanobis distance for all samples in the normal class is defined as: （11）； In equation (11): Indicates the first Weighted Mahalanobis distance of normal samples This indicates the first in the normal class. Standardized data for each sample; This represents the standard deviation of a normal sample. This represents the matrix transpose operation; ) is the feature weight matrix of the normal class samples; They represent the first The comprehensive weight of the first feature, the second feature The combined weight of each feature, No. The combined weight of each feature; （12）； In equation (12): This indicates the total number of ice-covered samples; （13）； The weighted Mahalanobis distance for all samples in the corresponding icing class is: （14）； In equation (14): Indicates the first Weighted Mahalanobis distance for each ice-covered sample; Indicates the first in the ice-covered category Standardized data for each sample.

5. The method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering as described in claim 4, characterized in that: The weighted dynamic Mahalanobis distance is calculated as follows: Assumption for A normal sample, in, express ; These represent the first normal sample and the second normal sample, respectively. ; Z These are standardized, normal sample data; Represents a vector augmented matrix; The augmented matrix of normal samples is expressed as follows , in, This represents the total augmented matrix at different lag times for normal samples; Let represent the augmented matrix at the first lag time of the normal sample, and represent the augmented matrix at the second lag time of the normal sample, respectively. Normal sample number An augmented matrix with a time lag; q Indicates the maximum lag order; Therefore, the weighted dynamic Mahalanobis distance (WDMD) for normal samples is: （15）； In equation (15): Indicates the first Weighted dynamic Mahalanobis distance of a normal sample This represents the standard deviation of a normal sample. （16）； In equation (16): Indicates the first in the normal class Standardized data for each sample; This represents the total number of normal samples; （17）； In equation (17): This represents the matrix transpose operation; Let the augmented matrix of the icing sample be... ; in: This represents the total augmented matrix at different lag times for the icing sample; Let represent the augmented matrix at the first lag time of the icing sample, and represent the augmented matrix at the second lag time of the icing sample, respectively. , frozen sample An augmented matrix with a time lag; The weighted dynamic Mahalanobis distance (WDMD) for the icing samples is calculated as follows: （18）； In equation (18): Indicates the first Weighted dynamic Mahalanobis distance of normal samples; （19）； In equation (19): This represents the standard deviation of a normal sample. This represents the total number of normal samples; Indicates the first in the normal class Standardized data for each sample; （20）； In equation (20): Indicates the first in the normal class Standardized data for each sample.

6. The method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering as described in claim 5, characterized in that: In step 3, the weighted dynamic Mahalanobis distance is mapped to the health status index of the wind turbine blades, including the following calculations: 3.1: Design of Icing Index for Wind Turbine Blades Based on Weighted Dynamic Mahalanobis Distance (WDMD) As shown in the following formula: （21）； In equation (21): Represented by natural constant e As the base, with and The product is an exponentiation operation; This indicates a regulatory index; Representing all samples The weighted dynamic Mahalanobis distance at time; because ,so ∈(0,1); （22）； In equation (22): This represents the confidence level under the condition that the wind turbine blades are healthy, normal, and not iced. The average of the normal samples WDMD ; （23）； In equation (23): The average of normal samples WDMD .

7. The method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering as described in claim 6, characterized in that: The WDMD of normal samples does not follow a normal distribution. In order to determine To obtain the value, the WDMD of normal samples must first be transformed into a normal distribution using appropriate transformations; the process of processing WDMD values ​​using the Box-Cox transformation is shown in the formula. The conversion formula is: （24）； In equation (24): The WDMD values ​​of normal samples after Box-Cox transformation. Represents a power transform parameter; It is estimated by the maximum likelihood method; First, construct the likelihood function. : （25）； In equation (25): yes Maximum likelihood estimation; Indicates a given When dealing with sample data, the likelihood value that the transformed sequence follows a normal distribution; This represents the WDMD value of a normal sample before the Box-Cox transformation; For all of equation (25) ,have: （26）； In equation (26): Indicates the first one sample W MD Take (λ) 1) Power; This represents the total number of normal samples; For each In Equation (26) Equation (27) can be expressed as: （27）； In equation (27): Normal samples after Box-Cox transformation The average value; （28）； In equation (28): This represents the total number of normal samples; Indicates the first in the normal class Standardized data for each sample; Index representing a normal sample; The following equation is derived: （29）； Before conversion It does not conform to a normal distribution. After Box-Cox transformation... The distribution characteristics satisfy the requirements of a normal distribution; using Box-Cox transformation and Criteria for determining thresholds As shown in equation (30), the confidence level of the wind turbine blades under healthy, normal, and non-icing conditions can then be determined. value; （30）； In equation (30): Normal samples after Box-Cox transformation The average value; For normal samples after Box-Cox transformation Standard deviation.

8. The method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering as described in claim 7, characterized in that: In step 4, a health status index threshold is calculated based on Gaussian mixture model (GMM) clustering, and then the health status index of the wind turbine blade icing state is quantitatively analyzed to obtain a quantitative classification of the severity of the icing state samples, specifically including: S4.1: Input data: Obtain the health status index and label reflecting the health level of the equipment from the sample: 0 indicates normal, 1 indicates icing; S4.2: Determine the optimal number of Gaussian clusters: Use the Bayesian Information Criterion (BIC) as the basis for model selection, as shown in the following equation (31): （31）； In equation (31): BIC(k) This represents the number of clusters determined based on the Bayesian information criterion. L It is the maximum value of the likelihood function corresponding to the model. K It is the number of clusters. n It is the sample size; As shown in equation (31), the greater the goodness of fit or the smaller the penalty value, the smaller the BIC value; therefore, the model with the smallest BIC value is selected as the optimal model; the number of clusters with the smallest BIC value is selected as the optimal number of clusters: （32）； In equation (32): , Both represent the optimal number of clusters.

9. The method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering as described in claim 8, characterized in that: Step 4 also includes GMM training: Using the optimal number of clusters Fit and train the GMM, then predict the cluster label and probability of each sample, and obtain the cluster label of each sample and the probability of each sample belonging to each Gaussian component: Initialize mean Covariance Weight In the iteration, each observation data is calculated. Belonging to the k The posterior probability of each probability distribution function is: （33）； In equation (33): Indicates the first i Sample Belonging to the k The posterior probability of each Gaussian component; For the first k The probability density function of a Gaussian distribution with Gaussian components; For the first k The mean of the distribution of the Gaussian components; For the first k Covariance of Gaussian components; For the first k The weights of each Gaussian component; Indicates the total number of clusters; （34）； In equation (34): Let represent the Gaussian kernel function, and represent the sample. With the k The center of each Gaussian component Similarity measure between them; Indicates the first i Observations of a sample; Indicates the first k The mean of the Gaussian components; Then, the following parameter values ​​are calculated: （35）； （36）； （37）； Get the mean of GMM Standard deviation and weight : （38）； In equation (38): They represent the first 1 The mean of the nth Gaussian components, the th 2 The mean of the Gaussian components, , No. k The mean of the Gaussian components; （39）； In equation (39): They represent the first 1 The covariance of the Gaussian components, the first 2 The covariance of the Gaussian components, , No. k Covariance of Gaussian components; （40）； In equation (40): They represent the first 1 The weight of the Gaussian component, the first 2 The weights of the Gaussian components, , No. k The weights of each Gaussian component; In equation (41): the Gaussian mixture model is used to obtain K One Gaussian component, parameters calculated by mean. Sort in ascending order, the corresponding standard deviation is The weight is ,and =1; In equation (42): for two adjacent Gaussian components and +1, we need to find a threshold. t This makes the weighted probability densities equal, where... Let be the probability density function of a Gaussian distribution. ; = In equation (43): Represents the Gaussian kernel function. Indicates the first k The mean of the Gaussian components, For the first k Covariance of Gaussian components; In equation (44): taking the natural logarithm of both sides of equation (43) and rearranging it into a quadratic equation; in a= （45）； In equation (45): a Indicates the first k The Gaussian component and the first k+ The difference in variance of one Gaussian component; For the first k The standard deviation of each Gaussian component; b= （46）； In equation (46): b Indicates the first k The Gaussian component and the first k+ The difference in the combined ratio of the mean to the variance of a single Gaussian component; Indicates the first k+ The mean of a single Gaussian component. Indicates the first k+ The standard deviation of one Gaussian component; c= （47）； In equation (47): c This indicates the influence of the combined weights and standard deviations of the two components on the location of the intersection point; Indicates the first k The weights of each Gaussian component; Indicates the first k+ The weight of a Gaussian component; Discriminant: （48）； In equation (48): The discriminant in the quadratic equation for threshold calculation in expression (44); if ≥0: （49）； （50）； In equations (49) and (50), This represents one of the solutions to the quadratic equation in expression (44). The expression (44) represents another solution to the quadratic equation; （51）； In equation (51), Indicates in and Choose a solution from among them.

10. The method for assessing the icing state of wind turbine blades based on dynamic weighted Mahalanobis distance and Gaussian mixture model clustering as described in claim 9, characterized in that: The calculated threshold is used as the quantitative standard for the health index of icing status, so as to realize the quantitative classification of the health index of icing status of wind turbine blades, thereby dividing the samples into different severity levels. According to the threshold, they are classified from low to high as ['severe icing', 'moderate icing', 'mild icing', 'minor icing']. （52）； Count the number of samples for each severity level: （53）； In formula (53): This represents the sample size for each severity level. This indicates the number of samples for each severity level. Indicates a statistical quantity.

Images