A wind turbine abnormality detection method and system based on SCADA periodic data

Abstract

This invention discloses a method and system for detecting anomalies in wind turbines based on SCADA periodic data. The method includes: acquiring operational data collected by the SCADA system of an actually operating wind turbine; constructing a normal dataset in addition to identifying abnormal data from SCADA alarms; the abnormal data collected by the SCADA system includes periodic, stable, and irregularly fluctuating data; initially labeling the three types of data in the abnormal data; using a neural network model to perform data pattern recognition on the acquired data to identify periodic data; performing feature engineering on the periodic data output to extract corresponding data features; enhancing the features using an isolation forest method to calculate anomaly scores; and using a machine learning model to classify the feature vectors formed within a time window to detect abnormal data within the periodic data. This method can improve the accuracy of anomaly detection in wind turbines.

Description

Technical Field

[0001] This invention relates to a method for detecting anomalies in wind turbine generators, specifically for detecting anomalies in wind turbine generator data. Background Technology

[0002] SCADA (Supervisory Control and Data Acquisition) is an indispensable core technology in modern wind power systems. By monitoring key operating parameters of wind turbine generators in real time, including wind speed, power generation, temperature, and vibration, it achieves comprehensive control over the equipment's operating status. The core functions of a SCADA system include fault detection and alarm, remote control and optimization, and data storage and analysis. Through long-term monitoring of wind turbine operating data, SCADA systems can support predictive maintenance, reducing downtime losses caused by equipment failures. In addition, SCADA systems provide environmental monitoring functions, dynamically adjusting wind turbine operating parameters by monitoring external environmental data such as wind speed and direction to optimize power generation efficiency and ensure safety. In large wind farms, SCADA systems can also achieve centralized management and control of multiple wind turbines, thereby improving the overall operating efficiency of the wind farm.

[0003] In actual wind power operation scenarios, anomalies can generally be categorized into three types: periodic, stationary, and irregularly fluctuating. Periodic anomalies are the most common, typically exhibiting time-related periodic fluctuations. For example, environmental parameters such as wind speed and direction usually fluctuate periodically with diurnal variations, while equipment operating parameters (such as power generation) may also show periodic changes under certain conditions. In contrast, stationary anomalies are usually characterized by some parameters deviating from normal values ​​for extended periods with relatively small fluctuations, such as power generation remaining consistently low over a certain period. Irregularly fluctuating anomalies are characterized by random fluctuations in parameters, which are difficult to predict and attribute, such as vibration data showing irregular and drastic changes at certain times. Among these three types of anomalies, periodic anomalies, due to their widespread prevalence and regularity, have become the focus of research. It is important to note that detecting periodic anomalies is not the same as anomaly detection itself; rather, it requires first distinguishing data with periodic characteristics from various data types and then designing detection methods based on these periodic characteristics. CN109086793A discloses an anomaly identification method for wind turbines. This method uses a series of data preprocessing and convolutional neural network methods to perform binary classification of wind turbine anomalies. However, while convolutional neural networks can extract local features well, they are not good at extracting temporal features. The use of a relatively complex network structure during online prediction is not conducive to real-time prediction. In addition, all anomalies are treated homogeneously, which fails to reflect the degree of anomaly and to assess the severity.

[0004] Meanwhile, existing studies generally suffer from insufficient feature extraction for wind turbine operating data, which exhibits significant periodicity. In actual production, wind turbine anomalies are closely related to the periodic changes in data. However, existing feature selection methods often employ simple statistical indicators, failing to fully consider the periodicity of the data, resulting in insufficient feature representation capabilities and ultimately affecting the accuracy of anomaly detection. Summary of the Invention

[0005] The technical problem to be solved by this invention is: how to identify periodic data patterns through neural networks, extract features from periodic abnormal data of wind turbines through feature engineering, automatically adjust the detection parameters of abnormal data of wind turbines, and improve the detection accuracy.

[0006] To address the aforementioned technical problems, this invention provides a method for detecting anomalies in wind turbines based on SCADA periodic data, comprising the following steps:

[0007] Step 1: Obtain operational data collected by the actual operating SCADA system of the wind turbine;

[0008] Step 2: Using all the collected operational data, construct a normal dataset by extracting abnormal data from SCADA alarms. The abnormal data collected by the SCADA system includes periodic, stable, and irregularly fluctuating data.

[0009] Step 3: Perform preliminary labeling of the three types of abnormal data: periodic, stationary, and irregularly fluctuating.

[0010] Step 4: Construct a neural network model and use the neural network model to perform data pattern recognition on the data obtained in Step 3 to identify periodic data.

[0011] Step 5: Perform feature engineering on the periodic data output from Step 4 to extract the corresponding data features;

[0012] Step 6: After performing multi-level feature engineering on the periodic data, the features are enhanced by calculating outlier scores using isolated forests.

[0013] Step 7: Build a machine learning model and use the machine learning model to classify the feature vectors formed within the time window for anomalies, and detect abnormal data in periodic data.

[0014] The aforementioned wind turbine anomaly detection method based on SCADA periodic data, in step 3, uses a phase division method integrating change point detection to initially label three types of data: periodic, stable, and irregularly fluctuating, including:

[0015] 31) A window-based change point detection method is used to discover potential change points; a bottom-up segmentation method is used to detect change points; and a pruned exact linear time algorithm is used to detect change points.

[0016] Three sets of change points are obtained by using the three algorithms mentioned above, and then a unified set of change points is obtained by using a change point fusion strategy based on hierarchical clustering.

[0017] Let the set of change points detected by the i-th method be . Where i∈{1,2,3}, n i Let be the number of change points detected by the i-th algorithm;

[0018] First, the sets of change points obtained from the three algorithms are merged, represented as follows:

[0019]

[0020] in To represent the total number of variable points after merging, the algorithm source identifier for each variable point is recorded as λ. r ∈{1,2,3}, r∈{1,2,...,N};

[0021] For each algorithm i, calculate the adjacent variable point spacing sequence D. i ={d i,j |d i,j =v i,j+1 -v i,j ,j∈[1,n i -1]};

[0022] Obtain the minimum variable point spacing

[0023] The time threshold TT is determined as TT = s × min i∈{1,2,3} (d i,min );

[0024] Where s∈[0.1,0.2] is the threshold coefficient, determined based on empirical data, v i,j This represents the j-th change point in algorithm i;

[0025] Construct the distance matrix D in this way p,q ∈R N×N D p,q =dist(v p ,v q );

[0026] In the hierarchical clustering process, any two variable points v are defined. p and v q The time distance between them is

[0027] dist(vp ,v q )=|v p -v q |

[0028] Hierarchical clustering is performed using the single-linkage criterion, and the inter-cluster distance is defined as:

[0029]

[0030] When D(C) i C j When )≤TT, merge cluster C i and cluster C j For cluster C formed after clustering k , containing m k A variable point

[0031] When m k When the value is ≥2, it indicates that at least two change points are generated within the cluster. In this case, the final change point location is determined by weighted averaging.

[0032]

[0033] Among them, w j For the weight, α i Let be the reliability coefficient of the i-th algorithm, and be a hyperparameter, both determined through multiple experiments. k c represents the final determined change point location in the k-th cluster. k,j Let j be the position of the variable point in the k-th cluster. To detect the variable point c k,j The reliability coefficient of the algorithm;

[0034] When m k When = 1, it indicates that only one change point detection algorithm has voted within the cluster, but it cannot determine whether the change point should be retained. In this case, by adjusting the time window [c k,1 -w,c k,1 Perform a statistical test using [+w], calculate the sample mean before and after the window, and calculate the pooled standard deviation, where w is a fixed value representing the number of windows.

[0035]

[0036] Where, σ pooled To combine standard deviation, x ind The value corresponding to the ind index, where ind is the coordinate index;

[0037] The statistical measure is used to assess the significance of the difference between the means before and after the window, reflecting the relative magnitude of the change in the means before and after the change point. The statistical measure is expressed as:

[0038]

[0039] D represents a standardized measure of mean change relative to data volatility;

[0040] Hypothesis testing is performed at a significance level α. If the probability value p, which measures the significance of the difference between the means, is less than α, then the variable point is retained, resulting in the final set of variable points P.

[0041] P = {p k |k∈[1,K]}

[0042] Where p k Let K represent the final coordinates of the variable points, and K be the final number of variable points, satisfying the condition...

[0043] min k∈[1,K-1] (p k+1 -p k )>TT

[0044] Calculate the interval ΔT between the points of change:

[0045] ΔT={Δt1,Δt2,...,Δt n-1}

[0046] Δt i =t i+1 -t i i = 1, 2, ..., n-1

[0047] t i Let represent the time corresponding to the i-th change point, and n represent the number of change points;

[0048] Calculate the interval variance and mean:

[0049]

[0050] Assume the ideal period of the data is T. ref The variance of the periodic data interval is The criteria for determining periodic data are:

[0051] And |μ ΔT -T ref |<ε μ ,

[0052] Where ε σ and ε μ are the tolerance thresholds for variance and mean, respectively, and is a hyperparameter.

[0053] The aforementioned wind turbine anomaly detection method based on SCADA periodic data, in step 4, constructs a CNN-BiLSTM hybrid neural network model and uses this model to perform data pattern recognition on the data obtained in step 3, identifying periodic data, including:

[0054] Step 41) Design a multi-layered CNN-BiLSTM hybrid neural network model, unfolding from input to output in the following order:

[0055] Input layer;

[0056] The first convolutional module includes convolutional layer one and pooling layer one;

[0057] The second convolutional module includes a second convolutional layer and a second pooling layer;

[0058] Bidirectional long short-term memory network layer;

[0059] Fully connected layer;

[0060] Output layer;

[0061] Step 42) The input layer receives normalized SCADA multidimensional time-series data, which is then processed by a one-dimensional CNN for feature extraction. During feature extraction, the CNN employs a sliding window mechanism, using convolution operations to extract local features.

[0062]

[0063] Where ki is the kernel size, w i Here, h represents the kernel weights, ReLU represents the activation function, b is the bias term used to adjust the output value of the convolutional layer; t The convolution output value at time step t is the feature generated after extracting local patterns from the input data;

[0064] Max pooling reduces feature dimensionality while preserving salient features;

[0065] The extracted feature sequence is input into a BiLSTM layer, including the forward pass. and reverse Two-way LSTM;

[0066]

[0067] x t This represents the feature vector input at time step t;

[0068] Finally, the classification result is output through a fully connected layer and a softmax function, using cross-entropy as the loss function.

[0069]

[0070] L is the total loss value, used to measure the difference between the predicted and the true label, nl is the total number of samples, and y i It is the true label of sample i. It is the probability distribution predicted by the model for sample i;

[0071] Step 43) The output of the CNN-BiLSTM model is a three-class classification result, which divides the data into periodic, stationary and irregular fluctuation types. For each data sequence, the model outputs the probability values ​​of periodic, stationary and irregular fluctuation types respectively.

[0072] The sum of all probability values ​​is 1, as shown in the following formula:

[0073] p periodic +p stationary +p irregular =1(1.16)

[0074] p periodic p represents the probability that the input belongs to a periodic category; stationary p represents the probability that the input belongs to a stationary category; irregular This represents the probability that the input belongs to the irregular fluctuation category;

[0075] The data type determination criterion is that the type with the highest output probability is the current type, and the corresponding formula is as follows:

[0076]

[0077] Where, p i For each probability type, p, s, and i represent periodic, stationary, and non-fluctuating types, respectively. The argmax function is used to obtain the parameters that maximize the input function value. The CNN-BiLSTM model directly outputs the type with the highest probability as the prediction result. For a data sequence identified as periodic, p... i =

[0078] p periodic Proceed to the next step of processing.

[0079] The aforementioned wind turbine anomaly detection method based on SCADA periodic data, in step 5, designs multi-level features, including simple thresholds and differentials. The simple threshold is a fixed threshold set for limit exceeding judgment, and the differential is the difference value calculated between adjacent times, days, and weeks.

[0080] Statistical methods include moving average and exponentially weighted moving average;

[0081] For data with obvious periodicity, temporal density analysis is used, with a one-week time window for density distribution analysis;

[0082] To address the abnormal patterns in multidimensional data, singular value decomposition is used for dimensionality reduction analysis to capture the underlying structure and main features of the data, and to identify the correlations and potential abnormal patterns between different dimensions.

[0083] Wavelet analysis is used to process anomalies at different frequency scales, and the week-on-week, day-on-day, and fluctuation ratio characteristics in the time domain are used to reflect periodic anomaly patterns.

[0084] The aforementioned wind turbine anomaly detection method based on SCADA periodic data includes, in step 6:

[0085] Step 61) Perform an isolation forest calculation on all features to obtain anomaly scores;

[0086] Taking week-on-week year-on-year comparison as an example, given time series data {X t For time t, the week-on-week year-on-year comparison is defined as the rate of change of the current value compared to the same period last week, WoW(t):

[0087]

[0088] x(t) represents the value of the week-on-week year-on-year characteristic at time t, and d represents the number of days;

[0089] Multiple window sizes W = {2, 3, 5, 7, 9} are selected to construct feature vectors;

[0090] For each window size w∈W, collect the week-on-week year-on-year sequence from the current time t back w weeks:

[0091]

[0092] f w (t) represents the week-on-week year-on-year feature vector based on window size w, containing the week-on-week year-on-year values ​​starting from the current time t and counting backwards for w weeks. All WoW(ti·7d) values ​​from i=0 to i=w-1 are arranged into a vector;

[0093] When w=3, the feature vector contains the year-on-year values ​​for the current week, the previous week, and the week before last:

[0094] f3(t)=[WoW(t),WoW(t-7d),WoW(t-14d)] (1.20)

[0095] By concatenating the feature vectors of all windows, we obtain the complete feature representation F(t):

[0096] F(t)=[f2(t),f3(t),f5(t),f7(t),f9(t)] (1.21)

[0097] For each time step in the time series, a similar feature vector is constructed to obtain the feature matrix X;

[0098] The outlier score of feature matrix X is calculated using the Isolation Forest algorithm:

[0099] For a sample point x, the anomaly score is defined as s(x,m):

[0100]

[0101] Where h(x) is the path length of the sample in the isolated tree, E(h(x)) is the average path length in all trees, and c(m) is the normalization factor:

[0102]

[0103] H(i) is the i-th harmonic number, H(i) = ln(i) + 0.5772156649, and m represents the sample size;

[0104] Finally, at each time t, an anomaly score s(t) is obtained. The anomaly score ranges between (-1, 1). The closer the anomaly score is to -1, the more abnormal it is. The closer it is to 1, the more normal it is.

[0105] Points with scores lower than a set value are identified as outliers;

[0106] For each feature window size w∈W, an anomaly score is obtained, for a total of 5 scores:

[0107] Step 62) Concatenate all features processed by the anomaly forest to form the anomaly score feature matrix:

[0108] S(t)=[s2(t),s3(t),s5(t),s7(t),s9(t)](1.24).

[0109] The aforementioned wind turbine anomaly detection method based on SCADA periodic data, in step 7, constructs an XGBoost model and uses the XGBoost model to classify anomalies in the feature vectors formed within the time window, detecting anomalous data in the periodic data, including:

[0110] Let the input be an n-dimensional feature matrix X∈R after feature engineering. m×n Where m represents the number of samples, n is the feature dimension, and each sample x i Contains n eigenvalues ​​[x i1 ,x i2 ,...,x in Meanwhile, the label vector is Y = [y1, y2, ..., y].m ] T , where y i ∈{0,1} indicates whether the i-th sample is an anomaly;

[0111] Standardize and preprocess the input features:

[0112]

[0113] Where μ j and σ j are the mean and standard deviation of the j-th feature, respectively;

[0114] The preprocessed data was divided into training and testing sets in an 8:2 ratio.

[0115]

[0116] X train Let X be the feature matrix of the training set. test Y is the feature matrix of the test set. train Y is the label vector of the training set. test For the test set label vector, y m Let m be the label of the m-th sample;

[0117] The XGBoost model constructs K decision trees iteratively, and the predicted value for each sample is:

[0118]

[0119] Where f k This represents the k-th decision tree;

[0120] In each iteration h, the model optimizes the following objective function:

[0121]

[0122] in, The gradient is first order; L (h) Let f be the loss after h iterations. k (x i ') represents the current k-th tree pair of sample x i The predicted value of ′ It is a loss function For the previous round of forecasts The second-order partial derivatives are used to measure the current predictive performance of the model;

[0123] Let γ be the second-order gradient, γ be the leaf node complexity penalty coefficient, λ be the L2 regularization coefficient, T be the number of leaf nodes, and w be the number of leaf nodes. j This is the predicted value for the j-th leaf node;

[0124] An early stopping strategy is used during model training. Training is stopped when the performance on the validation set does not improve within a set number of consecutive iterations.

[0125] For any test sample x i ', output its anomaly probability p i for

[0126]

[0127] The final category prediction is achieved by setting a threshold θ, as follows:

[0128]

[0129] The aforementioned wind turbine anomaly detection method based on SCADA periodic data, in step 4, constructs a ResNet neural network model and uses the ResNet neural network model for data pattern recognition, including:

[0130] The data sequence is segmented and normalized.

[0131] Label each type of time series sample with a category label: 0 for periodic types, and 1 for slowly changing types and irregularly fluctuating types.

[0132] We adopted a ResNet model structure improved by one-dimensional convolution. The input of the model is a segmented sample of time series with a shape of (batch size, 1, window size). Low-order features are extracted through an initial one-dimensional convolutional layer, followed by multiple residual modules to extract deep-level features. In each residual module, the data first goes through two layers of one-dimensional convolution and batch normalization. The input and output are directly added through skip connections to form residual features. The output of the residual module goes through a global average pooling layer to reduce the feature dimensionality. Finally, the probability of periodic and non-periodic categories is output through a fully connected layer.

[0133] The aforementioned wind turbine anomaly detection method based on SCADA periodic data, the forward propagation process of the ResNet model is represented as follows:

[0134] First, the input data is processed through an initial convolutional layer to generate features:

[0135] x1=ReLU(BatchNorm(Conv1D(x)))

[0136] Where ReLU represents the activation function, Conv1D represents the one-dimensional convolution function, and x1 is the initial feature map;

[0137] Then it passes through the residual module, and the output of the residual module is as follows:

[0138] x out =ReLU(BatchNorm(Conv1D(x1)))+x1

[0139] Where x out This represents the feature map after skip connections;

[0140] Finally, the output features of the residual module are fed into a global average pooling layer, and class probabilities are generated through a fully connected layer:

[0141]

[0142] GlobalAvgPool1D represents a one-dimensional global average pooling function, and Linear represents a linear transformation function.

[0143] During model training, the cross-entropy loss function is used as the optimization objective, and its formula is as follows:

[0144]

[0145] Among them, y i For real labels, The class probabilities predicted by the model;

[0146] After the model is trained, classification prediction is performed on newly added time series data. The new data is first preprocessed in the same way as the training data, including standardization and segmentation. Each segment of data is then input into the model sequentially to generate class prediction results. Let's assume the prediction result for a certain segment is p. i The final classification of the entire time series is determined by combining the prediction results of all segments through majority voting. The formula for majority voting is:

[0147]

[0148] Where c represents the category, M is the number of segments, and I is an indicator function used to count the number of predictions for each category. The category with the highest number of votes is considered the classification result. The final classification result of the time series is represented by p. i Let be the predicted category of the i-th time series segment by the model.

[0149] The aforementioned wind turbine anomaly detection method based on SCADA periodic data, in step 7, constructs a LightGBM model and uses the LightGBM model to classify anomalies in the feature vectors formed within the time window, including:

[0150] Let the input be an n-dimensional feature matrix X∈R after feature engineering. m×nWhere m represents the number of samples, n is the feature dimension, and each sample x i Contains n eigenvalues ​​[x i1 ,x i2 ,...,x in Meanwhile, the label vector is Y = [y1, y2, ..., y]. m ] T , where y i ∈{0,1} indicates whether the i-th sample is an anomaly;

[0151] LightGBM's optimization objective is expressed as:

[0152]

[0153] in The loss function is used to measure the model's predicted values. Compared with the true value y i The difference between them, Ω(f) is the regularization term used to control model complexity, defined as:

[0154]

[0155] Where T is the number of leaf nodes, w j γ is the weight of the j-th leaf node, and γ and λ are regularization coefficients. In binary classification tasks, the commonly used loss function is binary cross-entropy.

[0156]

[0157] in For sample x i The predicted probability;

[0158] LightGBM uses the gradient boosting method to gradually build a decision tree. First, it initializes the model's predictions. The mean of the target labels:

[0159]

[0160] In the h-th iteration, calculate the negative gradient and second gradient for each sample:

[0161]

[0162] Indicates the current predicted value Below, the first-order partial derivative of the loss function with respect to the predicted value, This represents the second-order partial derivative of the loss function with respect to the predicted value;

[0163] The split point of the new tree is determined by maximizing the gain, assuming that the current split point divides the dataset into two parts I. L and I R The gain formula is:

[0164]

[0165] in and Let be the sum of the first-order gradients of the left and right child nodes, respectively. and Given the second-order gradient sum, the feature with the largest gain and its split point are selected to complete the split by traversing all features and split points;

[0166] After H iterations, the model's final prediction is:

[0167]

[0168] in It is the Sigmoid function, used to map the model's output values ​​to the probability range [0,1]. The binary classification result is determined by setting a threshold τ. k (.) represents the predicted value of the k-th tree;

[0169]

[0170] An output of 1 indicates an error, and an output of 0 indicates no error.

[0171] A computer device / apparatus / system includes a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method described above.

[0172] The beneficial effects achieved by this invention are as follows: The method of this invention identifies periodic data patterns by combining convolutional neural networks and bidirectional long short-term memory networks, extracts features from periodic abnormal data through feature engineering, enhances periodic features by abstracting features using isolated forests, and builds a training model through XGBoost, which can automatically adjust the abnormal data detection parameters and improve the accuracy of abnormal data detection for wind turbines. Attached Figure Description

[0173] Figure 1 This is a flowchart of the wind turbine anomaly detection method based on SCADA periodic data in Embodiment 1 of the present invention;

[0174] Figure 2 This is a diagram of the CNN-BiLSTM classification model in Embodiment 1 of the invention. Detailed Implementation

[0175] The technical solution of the present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0176] Example 1

[0177] This embodiment provides a method for detecting anomalies in wind turbines based on SCADA periodic data, including the following steps:

[0178] Step 1: Obtain operational data collected by the actual operating SCADA system of the wind turbine;

[0179] Step 2: Using all the collected operational data, construct a normal dataset by extracting abnormal data from SCADA alarms. The abnormal data collected by the SCADA system includes periodic, stable, and irregularly fluctuating data.

[0180] Step 3: Initially label the three types of data—periodic, stationary, and irregularly fluctuating—using a phase-based change point detection ensemble. Then, address the change point detection problem from the perspectives of local dissimilarity, layer-by-layer aggregation, and global optimization using three methods: window-based change point detection, bottom-up segmentation strategy, and pruned exact linear-time method. These methods include:

[0181] A window-based change point detection method is used to discover potential change points. The method is as follows:

[0182] Two windows slide along the time series, respectively calculating the mean and variance of the data within the window, and measuring the difference between the mean and variance. When the difference between the two windows exceeds a preset threshold, the point of change is identified.

[0183] A bottom-up segmentation method is used to detect change points. The method is as follows:

[0184] A cohesive strategy is adopted, starting from the finest granularity, treating each point as an independent segment, and gradually merging adjacent sub-segments until the set stopping criteria are met. The merging standard is set based on the difference between the fitting error after merging and the error before merging, and unimportant variable points are gradually deleted to finally form the optimal segmentation.

[0185] The Pruned Exact Linear Time algorithm is used to detect change points. The method is as follows: Pruned Exact Linear Time is a dynamic programming algorithm that accurately locates change points by minimizing the cost function containing error and penalty terms. It also significantly reduces computational complexity through pruning rules. When the cost of a candidate point exceeds the constant of the current optimal solution, the candidate point and its subsequent points will be pruned and will no longer be considered as future change point candidates.

[0186] The three different change point detection methods described above are complementary. The window-based change point detection method is computationally simple, while the bottom-up segmentation method can supplement key changes that may be missed in window-based detection, making it suitable for handling complex stage partitioning problems. The global optimality of pruning with exact linear time compensates for the local optimization problems of the first two methods and is suitable for scenarios with large data scales. Integrating them can cover different data characteristics and application requirements, and by leveraging their respective advantages, the accuracy and robustness of stage partitioning can be improved.

[0187] Data was divided into three sets of change points using a seven-day period. Three different methods were employed to obtain these sets. Then, a hierarchical clustering-based change point fusion strategy was used to obtain a unified set of change points, including:

[0188] Let the set of change points detected by the i-th method be . Where i∈{1,2,3}, n i This represents the number of change points detected by this method;

[0189] First, the sets of change points obtained from the three algorithms are merged, represented as follows:

[0190]

[0191] in To represent the total number of variable points after merging, record the algorithm source identifier λ for each variable point. r r ∈{1,2,3}, r∈{1,...,N};

[0192] For each algorithm i, calculate the adjacent variable point spacing sequence D. i ={d i,j |d i,j =v i,j+1 -v i,j ,j∈[1,n i -1]}, obtain the minimum variable point spacing The time threshold TT is determined as TT = s × min i∈{1,2,3} (d i,min ), where s∈[0.1,0.2] is the threshold coefficient, determined based on empirical data, v i,j This represents the j-th change point in algorithm i;

[0193] Construct the distance matrix D in this way p,q ∈R N×N D p,q =dist(v p ,v q );

[0194] Next, in the hierarchical clustering process, define any two variable points v p and v qThe time distance between them is

[0195] dist(v p ,v q )=|v p -v q |

[0196] Hierarchical clustering is performed using the single-linkage criterion, and the inter-cluster distance is defined as:

[0197]

[0198] When D(C) i C j When )≤TT, merge cluster C i and cluster C j For cluster C formed after clustering k , containing m k A variable point According to m k The value is determined using a different strategy, which is based on the number of change points and simply reflects the source of the change point determination.

[0199] When m k When the value is ≥2, it indicates that at least two change points are generated within the cluster. In this case, a weighted average method is needed to determine the final change point location.

[0200]

[0201] Among them, w j For the weight, α i Let be the reliability coefficient of the i-th algorithm, and be a hyperparameter, both determined through multiple experiments. k c represents the final determined change point location in the k-th cluster. k,j Let j be the position of the variable point in the k-th cluster. To detect the variable point c k,j The reliability coefficient of the algorithm.

[0202] When m k When = 1, it indicates that only one change point detection algorithm has voted within the cluster, but it cannot determine whether the change point should be retained. In this case, by adjusting the time window [c k,1 -w,c k,1 Perform a statistical test using [+w], calculate the sample mean before and after the window, and calculate the pooled standard deviation, where w is a fixed value representing the number of windows.

[0203]

[0204] Where, σ pooled To combine standard deviation, x ind The value corresponding to the ind index, where ind is the coordinate index;

[0205] The statistical measure is used to assess the significance of the difference between the means before and after the window, reflecting the relative magnitude of the change in the means before and after the change point. The statistical measure is expressed as:

[0206]

[0207] D represents a standardized measure of mean change relative to data volatility;

[0208] Hypothesis testing is performed at a significance level of α = 0.05. If the probability value p, which measures the significance of the difference between the means, is less than α, the variable point is retained. After the above processing, the final set of variable points P is obtained:

[0209] P = {p k |k∈[1,K]}

[0210] Where p k Let K represent the final coordinates of the variable points, and K be the final number of variable points, satisfying the condition...

[0211] min k∈[1,K-1] (p k+1 -p k )>TT

[0212] Since there is a periodicity between the change points of periodic data, while there is no periodicity between the change points of slow-moving and irregularly fluctuating data, based on this characteristic, we first calculate the interval ΔT between the change points:

[0213] ΔT={Δt1,Δt2,...,Δt n-1}

[0214] Δt i =t i+1 -t i i = 1, 2, ..., n-1

[0215] t i Let represent the time corresponding to the i-th change point, and n represent the number of change points;

[0216] Calculate the interval variance and mean:

[0217]

[0218] Assume the ideal period of the data is T. ref The characteristic of periodic data is the variance of the interval. Smaller; mean of the interval μ ΔT Approaching the ideal period T ref Therefore, the criterion for determining periodic data is:

[0219] And |μ ΔT -Tref |<ε μ ,

[0220] Where ε σ and ε μ These are the tolerance thresholds for variance and mean, respectively. In practice, the ideal period T... ref The approximation between two peaks or troughs is used to simplify the data labeling problem into a binary classification problem. By using different strategies, it is determined whether the data is periodic, thus completing the labeling of the dataset. This process does not involve neural networks or machine learning, and a tolerance threshold needs to be specified, but it reduces the labeling complexity compared to manual methods.

[0221] The phase division method based on change point detection integration integrates multiple change point detection technologies that focus on different aspects, providing a solution for data phase division. It can obtain a set of change points and then automatically label the data according to the properties of periodic, slow, and irregular fluctuations.

[0222] Step 4: Construct a CNN-BiLSTM hybrid neural network model. Use this model to perform data pattern recognition on the data obtained in Step 3, identifying periodic data, including:

[0223] Step 41) Design a multi-layered neural network, unfolding it from input to output in the following order:

[0224] Input layer;

[0225] The first convolutional module includes convolutional layer one and pooling layer one;

[0226] The second convolutional module includes a second convolutional layer and a second pooling layer;

[0227] Bidirectional long short-term memory network layer;

[0228] Fully connected layer;

[0229] Output layer.

[0230] The input layer receives a three-dimensional tensor of shape (batch_size, time_steps, features), where:

[0231] X∈R batch_size×time_steps×features (1.6)

[0232] Where X is the data input to the neural network, which is a three-dimensional tensor, batch_size is the number of samples in a batch of data, time_steps is the time series length of each sample, features is the feature dimension of each time point, and R is the set of real numbers.

[0233] The feature extraction layer contains two convolutional modules, each consisting of a convolutional layer and a pooling layer;

[0234] The parameters of the first convolutional module are:

[0235]

[0236] The parameters of the second convolutional module are

[0237]

[0238] The sequence modeling layer consists of two stacked BiLSTM layers:

[0239] Sequence modeling as a functional description of bidirectional long short-term memory network layers

[0240]

[0241] The output tensor shape of BiLSTM is (batc_size, 2×units), where 2 represents bidirectional concatenation;

[0242] The classification layer uses two fully connected layers for final classification. The first fully connected layer uses the ReLU activation function, and the second fully connected layer uses the Softmax activation function, as shown below:

[0243]

[0244] The Softmax function ultimately transforms the output into a probability distribution.

[0245]

[0246] Where x is the input feature vector, representing the input data of the model; y i It is the category label, indicating one of the categories predicted by the model; P(y i |x) represents the model's predicted class y given input x. i The probability z represents the probability value of the output being of that category; i The model predicts the category y. i The original output value calculated is the output value before the activation function is applied; It is the result of exponentializing the output value of each category, representing the relative importance of that category; It is the sum of the index values ​​of all categories, used to normalize the index value of each category, ensuring that the sum of the probabilities of all categories is 1, z. j For the model, for category y j The original output value.

[0247] Add regularized Dropout between each layer to prevent overfitting.

[0248] The main advantage of using the CNN-BiLSTM model for anomaly detection is that it can comprehensively utilize the characteristics of convolutional neural networks and bidirectional long short-term memory networks to improve the accuracy and efficiency of detection. The CNN part is responsible for capturing local features in the data, while the BiLSTM part handles the global dependencies of the sequence data, enabling the model to consider both short-term patterns and long-term trends in the time series. This structure is particularly suitable for processing multidimensional time series data generated by SCADA systems because it can effectively identify and classify periodic, stationary, and irregularly fluctuating anomalies.

[0249] Step 42) The input layer receives normalized SCADA multidimensional time-series data, which is then processed by a one-dimensional CNN for feature extraction. During feature extraction, the CNN employs a sliding window mechanism, using convolution operations to extract local features.

[0250]

[0251] Where ki is the kernel size, w i Here, h represents the kernel weights, ReLU represents the activation function, b is the bias term used to adjust the output value of the convolutional layer; t The convolution output value at time step t is the feature generated after extracting local patterns from the input data.

[0252] Max pooling reduces feature dimensionality while preserving salient features;

[0253] The extracted feature sequence is input into a BiLSTM layer, including the forward pass. and reverse Two-way LSTM;

[0254]

[0255] x t Let t represent the feature vector input at time step t. The gating mechanism of LSTM can solve the gradient vanishing problem in long sequence training. The input gate, forget gate, and output gate control the input of new information, the forgetting of historical information, and the input of the current state, respectively.

[0256] BiLSTM's bidirectional structure can utilize contextual information from both before and after the sequence. Finally, it outputs the classification result through a fully connected layer and a softmax function, using cross-entropy as the loss function.

[0257]

[0258] L is the total loss value, used to measure the difference between the predicted and the true label, nl is the total number of samples, and y i It is the true label of sample i. It is the model's predicted probability distribution for sample i.

[0259] Step 43) The output of the CNN-BiLSTM model is a three-class classification result, which divides the data into periodic, stationary and irregular fluctuation types. For each data sequence, the model outputs the probability values ​​of periodic, stationary and irregular fluctuation types respectively.

[0260] Since this is a three-class classification problem, the sum of the probability values ​​is 1, as shown in the following formula:

[0261] p periodic +p stationary +p irregular =1(1.16)

[0262] p periodic p represents the probability that the input belongs to a periodic category; stationary p represents the probability that the input belongs to a stationary category; irregular This represents the probability that the input belongs to the irregular fluctuation category.

[0263] The data type determination criterion is that the type with the highest output probability is the current type, and the corresponding formula is as follows:

[0264]

[0265] Where, p i For each probability type, p, s, and i represent periodic, stationary, and non-fluctuating types, respectively. The argmax function is used to obtain the parameters that maximize the input function value. The CNN-BiLSTM model directly outputs the type with the highest probability as the prediction result. For a data sequence identified as periodic, p... i =

[0266] p periodic Proceed to the next step of processing.

[0267] Step 5: Perform feature engineering on the periodic data output from Step 4 to extract the corresponding data features.

[0268] To address the characteristics of SCADA system data, a multi-level feature design is implemented. This feature includes a simple threshold and a differential feature. The simple threshold is a fixed threshold used for limit-breaking judgment, while the differential feature calculates the difference between adjacent times, days, and weeks. This feature can represent unusual patterns such as trend changes and sudden fluctuations in the data, making it suitable for capturing abnormal fluctuations in different time periods.

[0269] Statistical methods used include moving average (MA) and exponentially weighted moving average (EWMA), where the window size for MA is 10 to 50 points and the smoothing coefficient α for EWMA is between 0.1 and 0.9.

[0270] EWMA is a weighted smoothing mechanism that helps to better capture recent trends, especially suitable for processes with rapid changes and cyclical fluctuations.

[0271] For data with obvious periodicity, temporal density analysis (TSD) is used to analyze the density distribution with a one-week time window. The time window setting is suitable for industrial process data and can help subsequent methods automatically adapt to changes in different time windows, thereby improving the sensitivity to periodic changes.

[0272] To address anomalous patterns in multidimensional data, singular value decomposition (SVD) is used for dimensionality reduction analysis. The number of rows ranges from 10 to 50, and the number of columns is 3, 5, or 7. This compresses the data into a low-dimensional representation, which is used to capture the underlying structure and main features of the data, identify the correlations between different dimensions and potential anomalous patterns, and reduce computational complexity.

[0273] Furthermore, wavelet analysis was used to process anomalies at different frequency scales, with window sizes of 3, 5, and 7 days, and time-domain features such as week-on-week comparisons, day-on-day comparisons, and fluctuation ratios were introduced. Wavelet analysis helps identify anomaly characteristics at different frequencies, effectively separating high-frequency noise from low-frequency trends, thus more accurately capturing the occurrence of anomalous signals. Time-domain features provide relative changes in data, helping to capture periodic anomalies and deviations relative to historical data, and helping to identify patterns that, while not violating a single threshold, are abnormal at a relative level.

[0274] Step 6: After performing multi-level feature engineering on the periodic data, in order to enhance the ability of features to represent outliers and avoid the significant impact of certain feature noise on subsequent classification models, the features are enhanced by calculating outlier scores using isolated forests, including:

[0275] Step 61) Perform an isolation forest calculation on all features to obtain anomaly scores;

[0276] Taking week-on-week year-on-year comparison as an example, given time series data {X t For time t, the week-on-week year-on-year comparison is defined as the rate of change of the current value compared to the same period last week, WoW(t):

[0277]

[0278] x(t) represents the value of the week-on-week year-on-year characteristic at time t, and d represents the number of days;

[0279] To capture the changing features at different time scales, multiple window sizes W = {2, 3, 5, 7, 9} are selected to construct feature vectors;

[0280] For each window size w∈W, collect the week-on-week year-on-year sequence from the current time t back w weeks:

[0281]

[0282] f w (t) represents the week-on-week year-on-year feature vector based on window size w, containing the week-on-week year-on-year values ​​starting from the current time t and counting backwards for w weeks. All WoW(ti·7d) values ​​from i=0 to i=w-1 are arranged into a vector.

[0283] When w=3, the feature vector contains the year-on-year values ​​for the current week, the previous week, and the week before last:

[0284] f3(t)=[WoW(t),WoW(t-7d),WoW(t-14d)] (1.20)

[0285] By concatenating the feature vectors of all windows, we obtain the complete feature representation F(t):

[0286] F(t)=[f2(t),f3(t),f5(t),f7(t),f9(t)] (1.21)

[0287] For each time step in the time series, a similar feature vector is constructed to obtain the feature matrix X.

[0288] The Isolation Forest algorithm is used to calculate the outlier score of the feature matrix X. The principle is that outliers are often easier to isolate.

[0289] For a sample point x, the anomaly score is defined as s(x,m):

[0290]

[0291] Where h(x) is the path length of the sample in the isolated tree, E(h(x)) is the average path length in all trees, and c(m) is the normalization factor:

[0292]

[0293] H(i) is the i-th harmonic number, H(i) = ln(i) + 0.5772156649; m represents the sample size;

[0294] Finally, at each time t, an anomaly score s(t) is obtained. The anomaly score ranges between (-1, 1). The closer the anomaly score is to -1, the more abnormal it is. The closer it is to 1, the more normal it is.

[0295] Points with scores less than -0.6 are usually considered outliers;

[0296] For each feature window size w∈W, an anomaly score is obtained, for a total of 5 scores:

[0297] The anomaly scores reflect the degree of anomaly at different time scales.

[0298] Step 62) Concatenate all features processed by the anomaly forest to form the anomaly score feature matrix:

[0299] S(t)=[s2(t),s3(t),s5(t),s7(t),s9(t)](1.24).

[0300] When the relationships between features are complex and non-linear, the original features may not be sufficient to effectively distinguish between normal and outlier points. Directly using these features for subsequent training cannot capture the features of anomalous samples well. Adding anomaly scores to an isolation forest can help the model distinguish important features, especially when the dataset contains a lot of noise. Anomaly scores can serve as an additional signal, indicating which data points may not conform to the normal pattern.

[0301] Step 7: Construct an XGBoost model and use it to classify anomalies in the feature vectors formed within the time window, detecting anomalous data in periodic data, including:

[0302] Let the input be an n-dimensional feature matrix X∈R after feature engineering. m×n Where m represents the number of samples, n is the feature dimension, and each sample x i Contains n eigenvalues ​​[x i1 ,x i2 ,...,x in Meanwhile, the label vector is Y = [y1, y2, ..., y]. m ] T , where y i ∈{0,1} indicates whether the i-th sample is an anomaly;

[0303] First, the input features are standardized and preprocessed:

[0304]

[0305] Where μ j and σ j are the mean and standard deviation of the j-th feature, respectively;

[0306] The preprocessed data was divided into training and testing sets in an 8:2 ratio.

[0307]

[0308]

[0309] X train Let X be the feature matrix of the training set. test Y is the feature matrix of the test set. train Y is the label vector of the training set. test For the test set label vector, y m Let m' be the label of the m-th sample;

[0310] The XGBoost model constructs K decision trees iteratively, and the predicted value for each sample is:

[0311]

[0312] Where f k This represents the k-th decision tree;

[0313] In each iteration h, the model optimizes the following objective function:

[0314]

[0315] in, The gradient is first order; L (h) Let f be the loss after h iterations. k (x' i ) represents the current k-th tree pair of sample x i The predicted value of ′ It is a loss function For the previous round of forecasts The second-order partial derivatives are used to measure the current predictive performance of the model;

[0316] Let γ be the second-order gradient, γ be the leaf node complexity penalty coefficient, λ be the L2 regularization coefficient, T be the number of leaf nodes, and w be the number of leaf nodes. j This is the predicted value for the j-th leaf node.

[0317] An early stopping strategy is used during model training. Training is stopped when the performance on the validation set does not improve in 50 consecutive iterations.

[0318] After each iteration, the current AUC (Area Under the Curve) and log loss value are output to monitor the training process.

[0319] For any test sample x i ', output its anomaly probability p i for

[0320]

[0321] The final category prediction is achieved by setting a threshold θ. In this embodiment, the threshold θ is set to 0.5, which is expressed as:

[0322]

[0323] Example 2

[0324] A method for detecting anomalies in wind turbines based on SCADA periodic data includes the following steps:

[0325] Step 1: Obtain operational data collected by the actual operating SCADA system of the wind turbine;

[0326] Step 2: Using all the collected operational data, construct a normal dataset by extracting abnormal data from SCADA alarms. The abnormal data collected by the SCADA system includes periodic, stable, and irregularly fluctuating types.

[0327] Step 3: Perform preliminary labeling of the three types of abnormal data: periodic, stationary, and irregularly fluctuating.

[0328] Step 4: Construct a neural network model and use the neural network model to perform data pattern recognition on the data obtained in Step 3 to identify periodic data.

[0329] Step 5: Perform feature engineering on the periodic data output from Step 4 to extract the corresponding data features;

[0330] Step 6: After performing multi-level feature engineering on the periodic data, the features are enhanced by calculating outlier scores using isolated forests.

[0331] Step 7: Build a machine learning model and use the machine learning model to classify the feature vectors formed within the time window for anomalies, and detect abnormal data in periodic data.

[0332] In step 4, a ResNet neural network model is constructed and used for data pattern recognition.

[0333] Residual Neural Network (ResNet) is a deep convolutional neural network architecture. The core innovation of ResNet lies in the introduction of residual blocks. By adding "skip connections" between network layers, the input signal is directly passed to subsequent layers, thereby alleviating the gradient vanishing problem in deep networks. ResNet can effectively support deep networks with hundreds or even thousands of layers.

[0334] Assume that each category of data—periodic, irregularly fluctuating, and slow-moving—contains 10 samples, each sample is a 7-day data sequence, and the sampling frequency is one data point per minute. Therefore, each sample contains 10080 data points.

[0335] To accommodate model training, the data sequence is first preprocessed. Since each sample contains many data points, directly inputting it into the network would be computationally expensive. Therefore, the time series is segmented, with each hour's data (60 points) serving as a window, resulting in 168 windows per sample. Each window is treated as an independent data sample. After segmentation, to eliminate differences in amplitude between different samples, a normalization method is used to normalize each segment. The normalization formula is:

[0336]

[0337] Where x is the original data, μ and σ represent the mean and standard deviation of the time series, respectively. The normalized data can ensure the numerical stability and training efficiency of the model, and x′ is the standardized data point value.

[0338] Label each type of time series sample with a category label: 0 for periodic types, and 1 for slowly changing types and irregularly fluctuating types.

[0339] To address the characteristics of time series data, a ResNet model structure with improved one-dimensional convolution is selected. The model's input consists of segmented samples of the time series, with a shape of (batch size, 1, window size). The model extracts low-order features through an initial one-dimensional convolutional layer, followed by multiple residual blocks to extract deeper features. In each residual block, the data first undergoes two layers of one-dimensional convolution and batch normalization. Then, the input and output are directly added together through a skip connection to form residual features. The output of the residual block passes through a global average pooling layer to reduce the feature dimensionality. Finally, a fully connected layer outputs the probabilities of periodic and non-periodic data, as step 3 can only distinguish between periodic and non-periodic data.

[0340] Specifically, the forward propagation process of the model can be represented as:

[0341] First, the input data is processed through an initial convolutional layer to generate features:

[0342] x1=ReLU(BatchNorm(Conv1D(x)))

[0343] Where ReLU represents the activation function, Conv1D represents the one-dimensional convolution function, and x1 is the initial feature map;

[0344] Then it passes through the residual module, and the output of the residual module is as follows:

[0345] x out =ReLU(BatchNorm(Conv1D(x1)))+x1,

[0346] Where x out This represents the feature map after skip connections;

[0347] Finally, the output features of the residual module are fed into a global average pooling layer, and class probabilities are generated through a fully connected layer:

[0348]

[0349] GlobalAvgPool1D represents a one-dimensional global average pooling function, and Linear represents a linear transformation function.

[0350] During model training, the cross-entropy loss function is used as the optimization objective, and its formula is as follows:

[0351]

[0352] Among them, y i For real labels, The model predicts the class probabilities. The Adam optimizer is used for training with an initial learning rate of 0.001. In each round of model training, data is input into the model in batches. In each batch, the gradient is calculated and the model parameters are updated until the loss function converges. L is the loss function.

[0353] After the model is trained, classification prediction is performed on newly added time series data. The new data is first preprocessed in the same way as the training data, including standardization and segmentation. Each segment of data is then input into the model sequentially to generate class prediction results. Let's assume the prediction result for a certain segment is p. i The final classification of the entire time series is determined by combining the prediction results of all segments through majority voting. The formula for majority voting is:

[0354]

[0355] Where c represents the category, M is the number of segments, and I is an indicator function used to count the number of predictions for each category. The category with the highest number of votes is considered the classification result. The final classification result of the time series is represented by p. iLet be the predicted category of the i-th time series segment by the model.

[0356] By segmenting time series data and combining it with the residual structure of ResNet, this approach efficiently captures both local and global features within the time series. Furthermore, the majority voting strategy enhances the model's ability to judge the entire time series. This method is suitable for large-scale time series classification tasks and exhibits high robustness and scalability.

[0357] Example 3

[0358] A method for detecting anomalies in wind turbines based on SCADA periodic data includes the following steps:

[0359] Step 1: Obtain operational data collected by the actual operating SCADA system of the wind turbine;

[0360] Step 2: Using all the collected operational data, construct a normal dataset by extracting abnormal data from SCADA alarms. The abnormal data collected by the SCADA system includes periodic, stable, and irregularly fluctuating types.

[0361] Step 3: Perform preliminary labeling of the three types of abnormal data: periodic, stationary, and irregularly fluctuating.

[0362] Step 4: Construct a neural network model and use the neural network model to perform data pattern recognition on the data obtained in Step 3 to identify periodic data.

[0363] Step 5: Perform feature engineering on the periodic data output from Step 4 to extract the corresponding data features;

[0364] Step 6: After performing multi-level feature engineering on the periodic data, the features are enhanced by calculating outlier scores using isolated forests.

[0365] Step 7: Build a machine learning model and use the machine learning model to classify the feature vectors formed within the time window for anomalies, and detect abnormal data in periodic data.

[0366] In step 7, a LightGBM model is constructed, and the LightGBM model is used to classify anomalies in the feature vectors formed within the time window;

[0367] LightGBM is a gradient boosting-based machine learning framework that optimizes the objective function by progressively building decision trees, thereby achieving high-precision predictions in classification tasks.

[0368] Let the input be an n-dimensional feature matrix X∈R after feature engineering. m×nWhere m represents the number of samples, n is the feature dimension, and each sample x i Contains n eigenvalues ​​[x i1 ,x i2 ,...,x in Meanwhile, the label vector is Y = [y1, y2, ..., y]. m ] T , where y i ∈{0,1} indicates whether the i-th sample is an anomaly;

[0369] The optimization goal of LightGBM can be expressed as:

[0370]

[0371] in The loss function is used to measure the model's predicted values. Compared with the true value y i The difference between them, Ω(f) is the regularization term used to control model complexity, defined as:

[0372]

[0373] Where T is the number of leaf nodes, w j γ is the weight of the j-th leaf node, and γ and λ are regularization coefficients. In binary classification tasks, the commonly used loss function is binary cross-entropy.

[0374]

[0375] in For sample x i The predicted probability;

[0376] LightGBM uses gradient boosting to gradually build a decision tree. First, it initializes the model's predictions. The mean of the target labels:

[0377]

[0378] In the h-th iteration, calculate the negative gradient (i.e., the first derivative of the loss function with respect to the predicted value) and the second gradient for each sample:

[0379]

[0380] Indicates the current predicted value Below, the first-order partial derivative of the loss function with respect to the predicted value, This represents the second-order partial derivative of the loss function with respect to the predicted value.

[0381] The split point of the new tree is determined by maximizing the gain, assuming that the current split point divides the dataset into two parts I. L and I R The gain formula is:

[0382]

[0383] in and Let be the sum of the first-order gradients of the left and right child nodes, respectively. and Given the second-order gradient sum, the feature with the largest gain and its split point are selected to complete the split by traversing all features and split points;

[0384] After H iterations, the model's final prediction is:

[0385]

[0386] in It is the Sigmoid function, used to map the model's output values ​​to the probability range [0,1]. The binary classification result is determined by setting a threshold τ. k (.) represents the predicted value of the k-th tree;

[0387]

[0388] An output of 1 indicates an error, and an output of 0 indicates no error.

[0389] A computer device / apparatus / system includes a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the above-described method.

[0390] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0391] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0392] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0393] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the technical principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A wind turbine anomaly detection method based on SCADA periodic data, characterized in that, Includes the following steps: Step 1: Obtain operational data collected by the actual operating SCADA system of the wind turbine generator; Step 2: Using all the collected operational data, construct a normal dataset by extracting abnormal data from SCADA alarms. The abnormal data collected by the SCADA system includes periodic, stable, and irregularly fluctuating data. Step 3: Perform preliminary labeling of the three types of abnormal data: periodic, stationary, and irregularly fluctuating. Step 4: Construct a neural network model and use the neural network model to perform data pattern recognition on the data obtained in Step 3 to identify periodic data. Step 5: Perform feature engineering on the periodic data output from Step 4 to extract the corresponding data features; Step 6: After performing multi-level feature engineering on the periodic data, the features are enhanced by calculating outlier scores using isolated forests. Step 7: Build a machine learning model and use the machine learning model to classify the feature vectors formed within the time window for anomalies, and detect abnormal data in periodic data. 2.The wind turbine anomaly detection method based on SCADA periodic data according to claim 1, wherein, In step 3, the stage division method of change point detection ensemble is used to initially label three types of data: periodic, stationary, and irregularly fluctuating, including: 31) A window-based change point detection method is used to discover potential change points; a bottom-up segmentation method is used to detect change points; and a pruned exact linear time algorithm is used to detect change points. Three sets of variable points are obtained using the three methods described above, and then a unified set of variable points is obtained by using a variable point fusion strategy based on hierarchical clustering. Set the first kind of method to detect the obtained variable point set as , , the number of variable points detected by the method; First, the sets of variable points obtained from the three methods are merged, represented as follows: ； in To determine the total number of variable points after merging, record the algorithm source identifier for each variable point. , ; For each algorithm Calculate the sequence of adjacent variable point distances. ; Obtain the minimum variable point spacing ; Time threshold Determined as ; in For threshold coefficient, Representation Algorithm The first in A turning point; Construct the distance matrix in this way , ; In the hierarchical clustering process, any two change points are defined. and The time distance between them is ； Hierarchical clustering is performed using the single-linkage criterion, and the inter-cluster distance is defined as: ； when At that time, merge clusters and cluster For clusters formed after clustering ,Include A variable point ; when If this occurs, it indicates that at least two change points are generated within the cluster. In this case, the final change point location is determined by weighted averaging. ； ； in, As weight, For the first The credibility coefficient of each algorithm. For the first The final determined location of the change point in each cluster For the first The th cluster The position of the variable point. To detect change points The reliability coefficient of the algorithm; when This indicates that only one change point detection algorithm has voted within the cluster, but it cannot determine whether the change point should be retained. In this case, the time window is adjusted accordingly. Perform statistical tests, calculate the sample mean before and after the window, and calculate the pooled standard deviation. The number of windows is a fixed value. ； ； in, To combine standard deviations, For the corresponding The value of the index For coordinate index; The statistical measure is used to assess the significance of the difference between the means before and after the window, reflecting the relative magnitude of the change in the means before and after the change point. The statistical measure is expressed as: ； A standardized measure of mean change relative to data volatility; At the significance level The following hypothesis test is performed to measure the probability value of the significance of the difference in means. Value less than If the variable is retained, the final set of variable points is obtained. : ； in Indicates the final coordinates of the variable point. The final number of variable points satisfies ； Calculate the interval between variable points : ； ； Indicates the first Each variable point corresponds to a specific time. Indicates the number of variable points; Calculate the interval variance and mean: ； ； Assuming the ideal period for the data is The variance of the periodic data interval is The criteria for determining periodic data are: ； in and These are the tolerance thresholds for variance and mean, respectively.

3. The wind turbine anomaly detection method based on SCADA periodic data according to claim 1, characterized in that, In step 4, a CNN-BiLSTM hybrid neural network model is constructed. This model is then used to perform data pattern recognition on the data obtained in step 3, identifying periodic data, including: Step 41) Design a multi-layered CNN-BiLSTM hybrid neural network model, unfolding from input to output in the following order: Input layer; The first convolutional module includes convolutional layer one and pooling layer one; The second convolutional module includes a second convolutional layer and a second pooling layer; Bidirectional long short-term memory network layer; Fully connected layer; Output layer; Step 42) The input layer receives normalized SCADA multidimensional time-series data, which is then processed by a one-dimensional CNN for feature extraction. During feature extraction, the CNN employs a sliding window mechanism, using convolution operations to extract local features. (1.12) in The kernel size is [size]. Here, ReLU represents the convolution kernel weights. It is a bias term used to adjust the value of the convolutional layer output; For time step The convolution output value is a feature generated after extracting local patterns from the input data; Max pooling reduces feature dimensionality while preserving salient features; The extracted feature sequences are input into a BiLSTM layer, including forward pass... and reverse Two-way LSTM; (1.13) Indicates time step The feature vector input at that time; Finally, the classification result is output through a fully connected layer and a softmax function, using cross-entropy as the loss function. (1.14) This is the total loss value, used to measure the difference between the predicted and the true labels. It is the total number of samples. It is a sample The true label, It is the model on the sample The predicted probability distribution; Step 43) The output of the CNN-BiLSTM model is a three-class classification result, dividing the data into periodic, stationary, and irregularly fluctuating types. For each data sequence, the model outputs the probability values ​​for periodic, stationary, and irregularly fluctuating types, respectively, with the sum of each probability value being 1. The formula is as follows: (1.16) This indicates the probability that the input belongs to a periodic category; This represents the probability that the input belongs to a stationary category; This represents the probability that the input belongs to the irregular fluctuation category; The data type determination criterion is that the type with the highest output probability is the current type, and the corresponding formula is as follows: (1.17) in, For each type of probability value, , , These represent periodic, stationary, and irregular fluctuation types, respectively. The function is used to obtain the parameters that maximize the value of the input function. The CNN-BiLSTM model directly outputs the type with the highest probability as the prediction result, for data sequences identified as periodic. = Proceed to the next step of processing.

4. The wind turbine anomaly detection method based on SCADA periodic data according to claim 1, characterized in that, In step 5, multi-level features are designed, including simple thresholds and differentials. The simple threshold is a fixed threshold set for limit violation judgment, and the differential is the difference value calculated between adjacent time, day and week respectively. Statistical methods include moving average and exponentially weighted moving average; For data with obvious periodicity, temporal density analysis is used, with a one-week time window for density distribution analysis; To address the abnormal patterns in multidimensional data, singular value decomposition is used for dimensionality reduction analysis to capture the underlying structure and main features of the data, and to identify the correlations and potential abnormal patterns between different dimensions. Wavelet analysis is used to process anomalies at different frequency scales, and the week-on-week, day-on-day, and fluctuation ratio characteristics in the time domain are used to reflect periodic anomaly patterns.

5. The wind turbine anomaly detection method based on SCADA periodic data according to claim 1, characterized in that, Step 6 includes: Step 61) Perform an isolation forest calculation on all features to obtain anomaly scores; For week-on-week comparisons, given time series data { }, for time Year-on-year comparison is defined as the rate of change of the current value relative to the value in the same period last week. : (1.18) Indicates the year-on-year characteristics of the week. The value of time, It represents the sky; Multiple window sizes of size 𝑊 = {2, 3, 5, 7, 9} are selected to construct feature vectors; For each window size ∈ 𝑊, collect data from the current time. Push forward Weekly year-on-year sequence: (1.19) Indicates based on window size The week-on-week eigenvectors contain the features from the current time. Start, push forward Weekly year-on-year value, From arrive All The values ​​are arranged sequentially into a vector; When 𝑤 = 3, the feature vector contains the year-over-year values ​​for the current week, the previous week, and the week before last: (1.20) By concatenating the feature vectors of all windows, a complete feature representation is obtained. : (1.21) For each time step in the time series, a similar feature vector is constructed to obtain the feature matrix 𝑋; The outlier score of the feature matrix 𝑋 is calculated using the Isolation Forest algorithm: For a sample point x, the anomaly score is defined as follows: : (1.22) in It is the path length of the sample in the isolated tree. It is the average path length across all trees. It is the normalization factor: (1.23) It is the first One harmonic number, = + 0.5772156649, Indicates the number of samples; Ultimately, every moment An abnormal score was obtained. The abnormal score ranges between (−1, 1). The closer the abnormal score is to −1, the more abnormal it is, and the closer it is to 1, the more normal it is. Points with scores lower than a set value are identified as outliers; For each feature window size ∈ Each of these will receive an abnormal score, for a total of 5 scores: Step 62) Concatenate all features processed by the anomaly forest to form an anomaly score feature matrix. (1.24).

6. The wind turbine anomaly detection method based on SCADA periodic data according to claim 1, characterized in that, In step 7, an XGBoost model is constructed, and the XGBoost model is used to classify anomalies in the feature vectors formed within the time window, detecting anomalous data in periodic data, including: Let the input be the result after feature engineering. 3D feature matrix ,in Indicates the number of samples. For the feature dimension, each sample Include eigenvalues Meanwhile, the label vector is ,in Indicates the first Is each sample an anomaly? Standardize and preprocess the input features: (1.25) in and The first The mean and standard deviation of the dimensional features; The preprocessed data was divided into training and testing sets in an 8:2 ratio. (1.26) (1.28) (1.29) (1.30) For the feature matrix of the training set, For the feature matrix of the test set, For the training set label vectors, For the test set label vector, For the first The labels of each sample; XGBoost models are built iteratively. Each of the decision trees has a predicted value for a given sample: (1.27) in Indicates the first A decision tree; In each iteration In this model, the objective function for optimization is as follows: (1.31) in, It is the first-order gradient; Represented as iteration Losses during the round Indicates the current number Tree samples Predicted values It is a loss function For the previous round of forecasts The second-order partial derivatives are used to measure the current predictive performance of the model; It is a second-order gradient. This represents the complexity penalty coefficient for leaf nodes. The L2 regularization coefficient is... The number of leaf nodes. For the first The predicted values ​​for each leaf node; An early stopping strategy is used during model training. Training is stopped when the performance on the validation set does not improve within a set number of consecutive iterations. For any test sample Output its anomaly probability for (1.32) By setting a threshold The final category prediction is represented as: (1.33)。 7. The wind turbine anomaly detection method based on SCADA periodic data according to claim 1, characterized in that, In step 4, a ResNet neural network model is constructed, and the ResNet neural network model is used for data pattern recognition, including: The data sequence is segmented and normalized. Label each type of time series sample with a category label: 0 for periodic samples, and 1 for stationary and irregularly fluctuating samples. The model uses a one-dimensional convolutional improved ResNet model structure. The input to the model is segmented samples of the time series, with the following shape: The system extracts low-level features through an initial one-dimensional convolutional layer, followed by multiple residual blocks to extract deeper features. In each residual block, the data first undergoes two layers of one-dimensional convolution and batch normalization, and then the input is directly added to the output through skip connections to form residual features. The output of the residual block passes through a global average pooling layer to reduce the feature dimension, and finally outputs the probabilities of the periodic and non-periodic categories through a fully connected layer.

8. The wind turbine anomaly detection method based on SCADA periodic data according to claim 7, characterized in that, The forward propagation process of the ResNet model is represented as follows: First, the input data is processed through an initial convolutional layer to generate features: ； in, Indicates activation function, Represents a one-dimensional convolution function. This is the initial feature mapping; Then it passes through the residual module, and the output of the residual module is as follows: ； in This represents the feature map after skip connections; Finally, the output features of the residual module are fed into a global average pooling layer, and class probabilities are generated through a fully connected layer: ； This represents a one-dimensional function for global average pooling. Represents a linear change function; During model training, the cross-entropy loss function is used as the optimization objective, and the formula is: ； in, For real labels, The class probabilities predicted by the model; After the model is trained, classification prediction is performed on newly added time series data. The new data is first preprocessed in the same way as the training data, including standardization and segmentation. Each segment of data is then input into the model sequentially to generate class prediction results. Let's assume the prediction result for a certain segment is... The final classification of the entire time series is determined by combining the prediction results of all segments through majority voting. The formula for majority voting is: ； in, Indicates category, For the number of segments, This is an indicator function used to count the number of predictions for each category. The category with the highest number of votes is considered the classification result. Represents the final classification result of the time series. For the model to the first Predicted categories for time series segments.

9. The wind turbine anomaly detection method based on SCADA periodic data according to claim 1, characterized in that, In step 7, a LightGBM model is constructed, and the LightGBM model is used to classify anomalies in the feature vectors formed within the time window, including: Let the input be the result after feature engineering. 3D feature matrix ,in Indicates the number of samples. For the feature dimension, each sample Include eigenvalues Meanwhile, the label vector is ,in Indicates the first Is each sample an anomaly? The optimization objective of LightGBM is expressed as: ； in The loss function is used to measure the model's predicted values. Compared with the true value The differences between them The regularization term, used to control model complexity, is defined as: ； in The number of leaf nodes. It is the first The weight of each leaf node, and These are all regularization coefficients. In binary classification tasks, the commonly used loss function is binary cross-entropy: ； in For the sample The predicted probability; LightGBM uses the gradient boosting method to gradually build a decision tree. First, it initializes the model's predictions. The mean of the target labels: ； In the In each iteration, the negative gradient and second gradient of each sample are calculated: ； Indicates the current predicted value Below, the first-order partial derivative of the loss function with respect to the predicted value, This represents the second-order partial derivative of the loss function with respect to the predicted value; The split point of the new tree is determined by maximizing the gain, assuming that the current split point divides the dataset into two parts. and The gain formula is: ； in and Let be the sum of the first-order gradients of the left and right child nodes, respectively. and Given the second-order gradient sum, the feature with the largest gain and its split point are selected to complete the split by traversing all features and split points; After completion After 10 iterations, the model's final prediction is: ； in It's the Sigmoid function, used to map the model's output values ​​to a probability range. The binary classification results are obtained by setting a threshold. To determine, It is the first The predicted value of each tree; ； An output of 1 indicates an error, and an output of 0 indicates no error.

10. A computer device / equipment / system, characterized in that, It includes a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method as claimed in any one of claims 1-9.

Images