A method and system for identifying icing regions of a fan blade

Abstract

The application relates to the technical field of icing area identification, and discloses a fan blade icing area identification method and system. The method comprises the following steps: obtaining a fan blade original image, performing light balance correction on the original image to obtain a corrected image; extracting a roughness matrix of the corrected image and completing area segmentation to obtain an initial icing boundary; extracting internal ice layer point cloud data based on the boundary, extracting geometric convex features to obtain an icing scale span; combining a preset scale span threshold value, performing differentiated thickness analysis on the ice layer point cloud to obtain ice layer thickness data; performing multi-dimensional feature classification on the thickness data to obtain an icing degree, combining historical icing data to calculate an icing evolution rate; extracting an evolution duration and matching an evolution type to finally obtain a blade icing state result. The method can realize high-precision identification of the icing area on the fan blade under complex working conditions.

Description

Technical Field

[0001] This invention relates to the field of icing area identification technology, and in particular to a method and system for identifying icing areas on wind turbine blades. Background Technology

[0002] Currently, wind power is a core component of the clean energy system. Blade icing in low-temperature and high-humidity environments can cause mass imbalance and aerodynamic performance degradation, directly leading to reduced power generation efficiency and even equipment damage. Accurate identification of icing areas and degrees of icing is crucial for the safe operation of wind turbines in winter. With the development of industrial vision technology, image-based non-contact detection has gradually become the mainstream approach for blade icing monitoring. However, existing solutions still struggle to simultaneously meet the dual requirements of boundary localization and thickness analysis, limiting overall recognition accuracy.

[0003] In existing technologies, blade icing detection often indirectly infers the overall icing state through vibration sensors or power curve analysis. Some industrial vision-based solutions simplify icing identification to image classification or target detection tasks. From the perspective of detection mechanism, indirect detection methods can only determine whether the entire machine is iced, but cannot locate specific icing areas and distribution patterns. Visual detection solutions mostly focus on a binary determination of whether or not icing has occurred, without incorporating the physical characteristics of continuous icing evolution. They cannot distinguish the cross-scale differences between thin frost adhesion and thick ice accumulation, and are also unable to cope with the visual similarity interference from changes in lighting, surface stains, and light icing, further exacerbating the deviations in boundary localization and thickness determination.

[0004] In summary, existing technologies suffer from insufficient accuracy in identifying icing regions under complex working conditions due to the lack of a collaborative mechanism for precise location of icing boundaries and differential analysis of ice thickness across scales. Summary of the Invention

[0005] This invention provides a method and system for identifying icing areas on wind turbine blades, so as to achieve accurate identification of icing areas on wind turbine blades under complex operating conditions.

[0006] In a first aspect, to solve the above-mentioned technical problems, the present invention provides a method for identifying icing areas on wind turbine blades, comprising: Acquire the original image of the wind turbine blades, and perform illumination equalization correction on the original image to obtain the corrected image; Texture features are extracted from the corrected image to obtain a roughness matrix. Region segmentation is then performed based on the roughness matrix to obtain the initial icing boundary. Based on the initial icing boundary, internal point cloud is extracted to obtain ice layer point cloud data. Based on the ice layer point cloud data, geometric protrusions are extracted to obtain the icing scale span. If the icing scale span is greater than a preset scale span threshold, the ice layer point cloud data is used to calculate the ice thickness to obtain the first ice layer thickness as the ice layer thickness data; if the icing scale span is less than or equal to the preset scale span threshold, the ice layer point cloud data is used to perform frost thickness inversion to obtain the second ice layer thickness as the ice layer thickness data. The ice thickness data is classified using multidimensional features to obtain the degree of icing. The difference between the degree of icing and the pre-acquired historical icing data is calculated to obtain the icing evolution rate. The duration of the icing evolution is extracted based on the icing evolution rate, and the evolution type is matched based on the duration of the icing evolution to obtain the final icing state.

[0007] Preferably, the step of acquiring the original image of the wind turbine blades and performing illumination equalization correction on the original image to obtain a corrected image includes: Acquire the original image of the wind turbine blade, extract the original pixel matrix of the original image, and generate a grayscale distribution histogram; Illumination distance analysis is performed based on the grayscale distribution histogram to obtain the mapping function; The target pixel matrix is ​​obtained by non-linearly stretching the original pixel matrix using the mapping function. The image is reconstructed based on the target pixel matrix to obtain the corrected image.

[0008] Preferably, the step of extracting texture features from the corrected image to obtain a roughness matrix, and then performing region segmentation based on the roughness matrix to obtain the initial icing boundary, includes: Pixel-to-pixel frequency statistics are performed on the corrected image to obtain the gray-level co-occurrence matrix; Feature extraction is performed based on the gray-level co-occurrence matrix to obtain gray-level image features, and spatial mapping is performed on the gray-level image features to obtain local roughness; The roughness matrix is ​​obtained by spatially arranging the local roughness. The roughness matrix is ​​used to predict icing, generating an icing probability matrix; The initial icing boundary is obtained by performing boundary segmentation based on the icing probability matrix.

[0009] Preferably, the step of extracting internal point clouds based on the initial icing boundary to obtain ice layer point cloud data, and extracting geometric protrusions based on the ice layer point cloud data to obtain the icing scale span, includes: The area within the initial icing boundary is extracted to obtain the ice layer region. The coordinates of the ice layer region are then transformed to obtain ice layer point cloud data. Curvature filtering is performed on the ice point cloud data to obtain a set of geometric protrusions; Calculate the vertical distance between the set of geometric protrusions and the preset reference plane to obtain the height dataset; The maximum height difference is calculated based on the height dataset to obtain the icing scale span.

[0010] Preferably, if the ice scale span is greater than a preset scale span threshold, then the ice layer point cloud data is used to calculate the ice thickness to obtain a first ice layer thickness as ice layer thickness data, including: If the ice scale span is greater than a preset scale span threshold, the ice point cloud data is cropped to obtain a thick ice point cloud subset. The surface normal vector is calculated based on the thick ice point cloud subset to obtain the normal vector set; The set of normal vectors is compared with the preset reference contour normal vectors, local extreme points are filtered out, and aggregated to obtain the ice layer protrusion region; Density analysis was performed on the ice protrusion region to obtain the spatial distribution density matrix; The volume integral of the ice protrusion region is performed based on the spatial distribution density matrix to obtain the first ice layer thickness as the ice layer thickness data.

[0011] Preferably, if the icing scale span is less than or equal to a preset scale span threshold, then the frost thickness is inverted on the ice point cloud data to obtain a second ice layer thickness as ice layer thickness data, including: If the ice scale span is less than or equal to the scale span threshold, feature extraction is performed on the ice layer point cloud data to obtain the point cloud roughness matrix; Gradient frequency domain transformation is performed on the point cloud roughness matrix to obtain the texture frequency matrix; Element-wise analysis of the texture frequency matrix yields the frost coverage rate; The refractive offset is obtained by refraction derivation based on the frost coverage rate. The refraction offset is used to perform refraction inversion calculations to obtain the second ice layer thickness as ice layer thickness data.

[0012] Preferably, the step of performing multi-dimensional feature classification on the ice thickness data to obtain the degree of icing, and calculating the difference between the degree of icing and pre-acquired historical icing data to obtain the icing evolution rate, includes: The thickness difference matrix is ​​obtained by calculating the pixel-by-pixel difference of the ice layer thickness data; Multiple features are extracted from the thickness difference matrix to obtain a fused feature vector; The fused feature vectors are mapped to a high-dimensional hyperplane to obtain the decision boundary distance; The degree of icing is obtained by classifying the data based on the decision boundary distance. Based on the degree of icing, obtain historical icing data for the corresponding collection period; The icing evolution rate is obtained by extrapolating the icing degree and the historical icing data.

[0013] Preferably, the step of extracting the evolution duration based on the icing evolution rate, performing evolution type matching based on the evolution duration, and obtaining the final icing state includes: Based on the freezing evolution rate and the preset rate threshold, time series analysis is performed to obtain the evolution duration. The physical evolution type is obtained by matching the evolution type based on the duration of the evolution. The final freezing state is obtained by mapping the state according to the physical evolution type.

[0014] In a second aspect, the present invention provides a wind turbine blade icing area identification system, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the method described above.

[0015] Thirdly, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the method described above.

[0016] Compared with the prior art, the present invention has the following beneficial effects: (1) By performing limited contrast adaptive histogram equalization correction on the original image of the leaf, the interference of local dark and overexposure caused by shooting distance and illumination angle is eliminated. Then, the contrast and energy features of surface texture are extracted based on the gray-level co-occurrence matrix and mapped to local roughness to generate a roughness matrix. This overcomes the problem that light icing and surface stains have similar visual appearance and low feature recognition under complex illumination, and improves the segmentation accuracy of icing area.

[0017] (2) By registering the internal region of the initial icing boundary with the standard three-dimensional model of the blade to generate ice layer point cloud data, extract the set of geometric protrusions and calculate the height extreme difference to obtain the icing scale span. Based on this, the thickness of thick ice is calculated by volume integration and the thickness of thin frost is calculated by refraction inversion, respectively, to achieve differentiated quantitative identification of icing morphology across scales.

[0018] (3) This invention calculates the icing evolution rate by using multi-cycle icing degree time series data, and combines the rate threshold with the evolution duration to match the physical evolution type, identify the dynamic process of rapid freezing and output the final icing state. This solves the problem that the existing technology only focuses on static icing results and cannot capture the process of rapid icing deterioration. It can provide a more valuable early warning basis for the safe operation of wind turbines in winter. Attached Figure Description

[0019] Figure 1 This is a schematic flowchart of a method for identifying icing areas on wind turbine blades provided in the first embodiment of the present invention. Detailed Implementation

[0020] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] Reference Figure 1 The first embodiment of the present invention provides a method for identifying icing areas on wind turbine blades, comprising the following steps: S11, acquire the original image of the wind turbine blades, perform illumination equalization correction on the original image, and obtain the corrected image; S12, extract texture features from the corrected image to obtain a roughness matrix, and perform region segmentation based on the roughness matrix to obtain the initial icing boundary; S13, extract the internal point cloud based on the initial icing boundary to obtain ice layer point cloud data, and extract geometric protrusions based on the ice layer point cloud data to obtain the icing scale span; S14. If the icing scale span is greater than a preset scale span threshold, the ice layer point cloud data is used to calculate the ice thickness to obtain the first ice layer thickness as the ice layer thickness data; if the icing scale span is less than or equal to the preset scale span threshold, the ice layer point cloud data is used to perform frost thickness inversion to obtain the second ice layer thickness as the ice layer thickness data. S15, perform multi-dimensional feature classification on the ice thickness data to obtain the degree of icing, and calculate the difference between the degree of icing and the pre-acquired historical icing data to obtain the icing evolution rate; S16, extract the evolution duration based on the freezing evolution rate, perform evolution type matching based on the evolution duration, and obtain the final freezing state.

[0022] In step S11, the original image of the wind turbine blades is acquired, and the original image is subjected to illumination equalization correction to obtain a corrected image, including: Acquire the original image of the wind turbine blade, extract the original pixel matrix of the original image, and generate a grayscale distribution histogram; Illumination distance analysis is performed based on the grayscale distribution histogram to obtain the mapping function; The target pixel matrix is ​​obtained by non-linearly stretching the original pixel matrix using the mapping function. The image is reconstructed based on the target pixel matrix to obtain the corrected image.

[0023] It should be noted that the original images of the wind turbine blades are acquired by industrial cameras deployed on the top of the nacelle or by visual acquisition equipment mounted on inspection drones. The acquisition angle covers the entire span of the blade, and the resolution of a single image is no less than 1920 x 1080 pixels. The image format is compatible with 8-bit single-channel grayscale images and three-channel color images. The acquisition cycle is adjustable from 5 to 30 minutes depending on the wind turbine's operating conditions. When extracting the original pixel matrix of the original image, the color image is first mapped to grayscale levels from 0 to 255 using a standard weighted grayscale conversion formula, generating a two-dimensional grayscale pixel matrix with the same spatial dimensions as the original image. Each element in the matrix corresponds to the grayscale value of a pixel, with the dimension being grayscale level. The total number of pixels corresponding to each grayscale level is counted, and the pixels are arranged in ascending order of grayscale value to generate a grayscale distribution histogram. The histogram can intuitively present the brightness and darkness distribution characteristics of the entire image, providing a statistical basis for subsequent light interference determination. The data processing accuracy of this step is controlled within ±1 gray level, and the processing time for a single image is no more than 20 milliseconds, which can meet the response requirements of real-time on-site inspection.

[0024] It is worth noting that the illumination distance analysis is based on the cumulative probability calculation of the gray-level distribution histogram. The percentage of pixels is accumulated at each gray level to obtain the cumulative probability sequence corresponding to levels 0 to 255. The focus is on the total cumulative probability value of the low gray-level range (0 to 50), which corresponds to the darker areas of the image. Its proportion directly reflects the degree of insufficient illumination caused by excessive shooting distance or backlighting. The distribution judgment condition is set as a cumulative probability greater than 80% in the low gray-level range. This threshold aims to maximize the accuracy of icing recognition after illumination correction. ROC curve analysis was performed on 120 sets of leaf images covering multiple illumination scenarios. The cumulative probability corresponding to the maximum Youden's index was taken as 80% as the final threshold, which can accurately identify more than 92% of images with insufficient illumination. When the judgment condition is met, it indicates that the overall image is dark, and a mapping function related to the degree of illumination distance interference needs to be constructed. The function is based on contrast-limited adaptive histogram equalization, and its core expression is:

[0025] in The original pixel grayscale value. and These are the gray levels of the dark area boundaries and bright area boundaries dynamically extracted based on the cumulative probability of low gray levels. Specifically, the gray level corresponding to when the cumulative probability of the low gray level range reaches 80% is taken as... This value represents the maximum grayscale value of a pixel that is darker due to lighting distance interference; the lower the value, the more severe the lighting distance interference. The grayscale level corresponding to an overall cumulative probability of 95% is taken as... This represents the boundary of the bright areas in the image that still retain effective texture information. The mapping function linearly stretches pixels originally concentrated in a narrow low-grayscale range to the full grayscale range, significantly improving the grayscale values ​​of pixels in dark areas and specifically correcting localized darkening caused by illumination distance interference. The function also introduces a contrast limiting threshold of 2.0, performing local mapping in 8x8 pixel blocks to avoid excessive amplification of background noise, thus maintaining a natural transition of image details while removing illumination distance interference. This threshold setting is based on the common engineering specifications for limiting contrast adaptive histogram equalization in the field of industrial visual inspection, and has been calibrated through testing on 120 sets of leaf images. This parameter combination can improve dark details while keeping the noise amplification factor within 1.2 times, balancing enhancement effect and noise suppression. To avoid excessive amplification of background noise, the mapping function can expand the original narrow grayscale range to the full grayscale range, achieving dynamic adjustment of the brightness weights of different areas.

[0026] It's worth further explaining that the nonlinear stretching operation takes the original pixel matrix as input and substitutes it pixel by pixel into the constructed mapping function for grayscale conversion. Shadow areas with low grayscale values ​​in the original pixel matrix have their grayscale values ​​significantly increased after mapping, while pixels with grayscale values ​​in the normal range maintain a relatively smooth transition, avoiding light-dark banding. This nonlinear stretching process can expand the image's grayscale dynamic range by approximately 2.5 times, effectively highlighting the surface details of leaves hidden in low-light areas. Contrast limit verification is performed simultaneously during the stretching process. If the grayscale stretching amplitude in a local block is too large, causing the contrast to exceed the limit, the mapping coefficient for that area is attenuated to ensure that noise does not amplify synchronously with brightness. After the entire matrix is ​​stretched, a target pixel matrix with the exact same size as the original pixel matrix is ​​obtained. The calculation process for a single 1080P image can be completed within 15 milliseconds.

[0027] Subsequently, image reconstruction is performed based on the target pixel matrix. Following the spatial coordinate order of the original image, the grayscale values ​​within the matrix are mapped point-by-point to pixels, generating a corrected image in single-channel grayscale format. The reconstruction process maintains no offset in pixel spatial position, ensuring the image resolution is completely consistent with the original image. Sub-pixel-level coordinate alignment ensures that subsequent feature extraction accurately corresponds to the original acquisition perspective, with registration errors controlled within 0.5 pixels. The corrected image effectively eliminates interference from changes in illumination angle, shooting distance, and other factors such as localized underexposure and overexposure due to reflections. Fine texture features on the leaf surface, originally hidden in shadows or strong light, are clearly revealed. Verified with 120 sets of field samples, this correction step reduces the boundary positioning error of subsequent icing area segmentation by approximately 18% and increases image information entropy by an average of 32%. After the corrected image is generated, it is stored in a local image cache unit, serving as the sole input data source for subsequent texture feature extraction steps.

[0028] In step S12, texture features are extracted from the corrected image to obtain a roughness matrix. Region segmentation is then performed based on the roughness matrix to obtain the initial icing boundary, including: Pixel-to-pixel frequency statistics are performed on the corrected image to obtain the gray-level co-occurrence matrix; Feature extraction is performed based on the gray-level co-occurrence matrix to obtain gray-level image features, and spatial mapping is performed on the gray-level image features to obtain local roughness; The roughness matrix is ​​obtained by spatially arranging the local roughness. The roughness matrix is ​​used to predict icing, generating an icing probability matrix; The initial icing boundary is obtained by performing boundary segmentation based on the icing probability matrix.

[0029] It should be noted that pixel pair frequency statistics are performed based on a local sliding window of the calibrated image. The calibrated image is divided into multiple sliding calculation windows, with a window size of 11×11 pixels. This size is statistically calibrated using typical texture scales on the blade surface. Aiming for the optimal balance between texture classification accuracy and single-image computation time, five window sizes from 5×5 to 15×15 were tested. The 11×11 window achieved the highest overall score, balancing local texture capture capability with regional statistical stability. Within each window, a pixel step size of 1 and a statistical direction of 0 degrees are set, counting the frequency of occurrence of adjacent pixel pairs in the horizontal direction. All pixel combinations within the window are traversed, and the occurrence frequency of each grayscale value pair is counted. A two-dimensional grayscale co-occurrence matrix is ​​generated by arranging the grayscale levels, with 256 rows and columns, corresponding to the total number of grayscale levels. All elements in the grayscale co-occurrence matrix are divided by the total frequency of pixel pairs within the window to obtain a normalized co-occurrence matrix. The elements in this matrix range from 0 to 1, representing the joint occurrence probability of corresponding grayscale pixel pairs. The window sliding step size is set to 6 pixels. After all windows in the entire image are statistically analyzed, a set of gray-level co-occurrence matrices of the corresponding size is obtained.

[0030] It is worth noting that grayscale image features include two key statistical quantities: contrast and energy value. Contrast is calculated by weighted summation of the elements and grayscale differences within the grayscale co-occurrence matrix. A higher value indicates more drastic local grayscale changes and more pronounced surface texture undulations. Energy value is calculated by the sum of squares of matrix elements. A higher value indicates a more uniform grayscale distribution and a smoother surface. The preset grayscale level thresholds are a contrast of no less than 20 and an energy value of no more than 0.3. This threshold set is based on the texture statistics of 200 normal, ice-free leaves. The average contrast of normal, smooth leaves is approximately 12, and the average energy is approximately 0.45. An energy value lowered to 0.3 is used as the judgment boundary for deviations from the normal state. This threshold can effectively filter texture regions that deviate from the normal state. A spatial mapping operation is performed on grayscale image features that meet the threshold conditions to construct a mapping rule from contrast and energy combinations to roughness levels. The feature combinations are mapped to integer roughness levels from 1 to 5, with higher values ​​indicating a rougher surface. The roughness level output by a single window is the local roughness of the corresponding area, and the recognition accuracy of the thin frost area can reach more than 85%.

[0031] It is worth further explaining that the spatial arrangement of local roughness strictly follows the coordinate system of the original image. It traverses all local calculation windows within the corrected image, and fills the corresponding local roughness values ​​output by each window according to the row and column positions of the window center in the original image. Image edge regions are processed using a mirror filling method to avoid edge distortion, generating a two-dimensional roughness matrix proportional to the spatial resolution of the corrected image. The spatial position of each element in the matrix corresponds one-to-one with the corresponding region in the original image, without coordinate offset or scale scaling, ensuring that the texture distribution accurately matches the actual area of ​​the blade, with a registration error of less than one window unit. The values ​​in the roughness matrix are uniformly distributed in discrete levels from 1 to 5, intuitively presenting the spatial distribution pattern of the roughness of the blade surface texture. The roughness matrix preserves the spatial topological information of the image while weakening the illumination interference of the original grayscale. The feature discrimination between the icing area and the normal area is improved by about 30% compared to the original grayscale image, providing more discriminative input features for subsequent icing area segmentation.

[0032] Subsequently, icing prediction is performed using a pre-trained convolutional neural network. The network architecture employs a simplified residual network structure, comprising four stages. Each stage consists of two stacked residual modules, each containing two 3×3 convolutional layers (using the ReLU activation function) and an identity skip connection. Stages are downsampled via max-pooling layers with a stride of 2. The decoding path is symmetrically configured with four stages. Each stage first performs upsampling using a 2x bilinear interpolation, then skips the feature map with the corresponding stage of the encoding path, followed by two 3×3 convolutional layers to restore spatial resolution. The final output layer is a 1×1 convolutional layer with a sigmoid activation function, progressively restoring the low-resolution feature map to a 256×256 pixel-level probability output consistent with the input roughness matrix size. The network consists of a total of 25 convolutional layers and 4 max-pooling layers. The pre-training dataset contains 12,000 annotated images of wind turbine blade icing, covering various operating conditions with different lighting and icing levels. The annotations were cross-validated by three wind power operation and maintenance experts, achieving pixel-level accuracy. During the training phase, the Adam optimizer was used with an initial learning rate of 0.001 and a batch size of 16, for a total of 50 training epochs. The convolutional neural network was trained using pixel-level binary cross-entropy loss, and the segmentation effect of the icing region was evaluated based on the intersection-union ratio and Dice coefficient of the validation set. The network parameters corresponding to the epoch with the best segmentation effect on the validation set were saved, and the classification accuracy on the validation set reached 92.3%. The roughness matrix was uniformly adjusted to a size of 256 by 256 and used as the network input. Deep semantic features were extracted layer by layer through multiple convolution and pooling operations. The output layer obtained the unnormalized value of the icing category corresponding to each pixel. Then through the Sigmoid function The probability values ​​mapped to 0 to 1 are multiplied by 100 to generate an icing probability matrix. The elements in the matrix range from 0 to 100, with higher values ​​indicating a greater likelihood that the corresponding pixel is an icing area.

[0033] It should be noted that boundary segmentation first performs binarization processing based on the icing probability matrix, with a preset segmentation threshold of 85. This threshold is calibrated based on a PR curve plotted from 3000 pixel-level labeled validation samples, using the probability value of 85 corresponding to the maximum F1 score of 0.89 as the segmentation threshold, which can simultaneously ensure both precision and recall. All pixels within the probability matrix are traversed; if a pixel's probability value is greater than the preset segmentation threshold, it is marked as an icing region and assigned a value of 1; otherwise, it is marked as a non-icing region and assigned a value of 0, thus generating an initial segmentation mask in binary format. A 3×3 Sobel edge detection operator is used to perform convolution operations on the initial segmentation mask, calculating the grayscale gradient magnitudes in the horizontal and vertical directions respectively. The gradient magnitudes are normalized to the 0-1 range and used as the edge weights for the corresponding pixels; a higher weight value indicates a higher probability that the pixel is at the boundary between icing and non-icing. Guided by edge weights, this method refines the segmentation boundary, separating two types of regions by minimizing the boundary segmentation cost. Specifically, the edge weights are used as grayscale values ​​in the terrain height map, and the icing probability mask is used as the initial marker. A watershed transformation is applied to expand the edge weight gradient map until the two types of markers meet, forming a closed segmentation line. This segmentation line is the initial icing boundary, represented by a closed polygon with vertices representing a series of consecutive pixel coordinates and a boundary line width of 1 pixel. The boundary coordinates can be directly mapped back to the original image space, achieving pixel-level (±1 pixel) boundary localization accuracy. Testing shows that this segmentation method achieves a boundary overlap of 0.89, significantly better than the 0.72 of direct grayscale segmentation, accurately defining the contour range of the icing region and providing a 2D contour basis for subsequent 3D point cloud extraction.

[0034] In step S13, internal point cloud extraction is performed based on the initial icing boundary to obtain ice layer point cloud data. Geometric protrusion extraction is then performed based on the ice layer point cloud data to obtain the icing scale span, including: The area within the initial icing boundary is extracted to obtain the ice layer region. The coordinates of the ice layer region are then transformed to obtain ice layer point cloud data. Curvature filtering is performed on the ice point cloud data to obtain a set of geometric protrusions; Calculate the vertical distance between the set of geometric protrusions and the preset reference plane to obtain the height dataset; The maximum height difference is calculated based on the height dataset to obtain the icing scale span.

[0035] It should be noted that the internal region of the initial icing boundary is extracted based on the pixel coordinate system of the blade image. A polygon filling method is used to retain all pixel regions within the boundary contour, while removing non-iced pixels on the outside. During the acquisition of the original image, depth data corresponding to each pixel is recorded simultaneously. This depth data is acquired by a laser ranging module or binocular stereo vision unit mounted on the inspection equipment, and is aligned pixel-by-pixel with the RGB image, with a depth measurement error of no more than ±2 mm. Spatial mapping uses the standard 3D model of the blade as the geometric reference. This 3D model is provided by the wind turbine manufacturer and is in the STEP general 3D format. A coordinate system is established with the center of the flange at the blade root as the origin, the spanwise direction as the X-axis, the chordwise direction as the Y-axis, and the thickness direction as the Z-axis. The model surface accuracy is ±1 mm. The pixels inside the boundary are downsampled at equal intervals, with the center pixel taken from every 3×3 pixel block. For each downsampled pixel, its 3D position in the camera coordinate system is determined based on its corresponding measured depth value. This position is then transformed to the blade model coordinate system using a pre-calibrated camera extrinsic matrix, yielding the 3D discrete point coordinates of the actual ice layer's outer surface. The set of all discrete points constitutes the ice layer point cloud data. The standard 3D blade model is only used to provide a coordinate registration benchmark and a reference surface in the ice-free state; it does not replace the actual ice layer depth. The depth data is primarily used to determine the 3D spatial location, geometric protrusions, and ice scale span of the icing area. For thin frost areas with a thickness less than the depth measurement error, the frost thickness is not directly determined based on the point cloud elevation difference. Instead, the spatial range of the thin frost area is determined using the point cloud data, and the second ice layer thickness is inverted in step S14 based on the image refraction offset of the same area. Point cloud data contains the X, Y, and Z spatial coordinates of each point, measured in millimeters. The point cloud density is approximately 2 points per square millimeter. At this density, a 5mm grid can contain about 50 points, sufficient to support the eigenvalue decomposition of the neighborhood covariance matrix and stably calculate local curvature. The number of point clouds in a single icing region is approximately 30,000 to 60,000, which can completely represent the three-dimensional topological structure of the icing region, providing a three-dimensional data source for subsequent geometric protrusion extraction.

[0036] Curvature selection first calculates the local curvature features of each data point in the ice layer point cloud data. A neighborhood set is formed by selecting the 20 nearest neighbors of each point, with a point cloud density of 2 points per square millimeter. The 20 neighborhood points correspond to a local computational region of approximately 10 square millimeters. This neighborhood size can stably fit the local surface normal vector and curvature while suppressing interference from single-point noise. The local curvature value is calculated based on the surface normal variation amplitude of the neighborhood points. A covariance matrix is ​​constructed for each point and its 20 nearest neighbors. The eigenvalues ​​of the covariance matrix are calculated, and the ratio of the smallest eigenvalue to the sum of the eigenvalues ​​is the local curvature value of that point, with a dimension of 1. The preset dimensionless curvature threshold is 0.25. This threshold calibration is based on the statistical calibration of the curvature of labeled ice protrusions and normal leaf surfaces in the ice point cloud. The ratio of covariance eigenvalues ​​is used as the curvature measure. The larger the ratio, the more drastic the surface change. This can distinguish between the normal leaf surface and ice protrusions while excluding gentle surface changes caused by slight frost adhesion, ensuring that only significant geometric protrusions such as ice nodules and ice ridges are included in the set. All point cloud data are traversed, and data points with curvature values ​​greater than the preset curvature threshold are extracted and added to the geometric protrusion set. Each element in the set contains the three-dimensional coordinates of that point. The corresponding curvature values ​​are used for subsequent analysis of protrusion height. This set corresponds to irregular protrusion structures such as ice nodules and ice ridges on the ice surface. Verified by field measurements of ice nodule samples, the number of geometric protrusion sets under this threshold accounts for 8% to 20% of the total ice point cloud, which is basically consistent with the proportion of protrusion area observed by physical observation.

[0037] It is worth noting that the normal un-iced surface region adjacent to the initial icing boundary in the standard 3D model of the blade was selected as the reference plane selection area, and the plane equation was directly constructed by extracting points within this selection area from the point cloud data. The parameters of the reference plane are obtained by fitting using the least squares method, and this is used as the preset reference plane. The reference plane is only used for the local neighborhood corresponding to the current icing area, and the fitting result is constrained by the local tangent plane direction of the standard three-dimensional model of the blade at that location to avoid the blade's own curvature and torsion being included in the ice protrusion height. The perpendicular distance from all data points within the geometric protrusion set to the reference plane is calculated point by point. For any point... Its plane The vertical distance is

[0038] All calculated distance values ​​are arranged according to their spatial location in the point cloud, forming a height dataset. The data unit is millimeters, and the numerical range within the dataset corresponds to the height distribution of ice layer protrusions. Verified using multiple sets of standard gauge blocks, the distance calculation error can be controlled within ±1 millimeter, meeting the accuracy requirements for ice layer thickness analysis.

[0039] It should be further explained that the icing scale span is calculated based on the height dataset. All distance values ​​within the dataset are sorted, and the maximum and minimum height values ​​are selected. The maximum height corresponds to the highest point of the ice layer protrusion, and the minimum height corresponds to the baseline bottom position of the protrusion area, i.e., close to the blade body. Subtracting the minimum height from the maximum height gives the maximum height difference, which is the icing scale span. The icing scale span is measured in millimeters and comprehensively characterizes the degree of three-dimensional morphological difference of the icing area in the direction perpendicular to the blade surface. A larger span indicates more dramatic morphological fluctuations in the ice layer and a higher probability of thick ice accumulation. For example, when the extreme height values ​​are 45 mm and 12 mm, the icing scale span is 33 mm. This parameter is a key criterion for subsequent differential thickness analysis, enabling automatic separation of thick ice and thin frost scenarios, with a relative error of less than 5% compared to manual measurements.

[0040] In step S14, if the ice scale span is greater than a preset scale span threshold, the ice layer point cloud data is used to calculate the ice thickness to obtain the first ice layer thickness as ice layer thickness data, including: If the ice scale span is greater than a preset scale span threshold, the ice point cloud data is cropped to obtain a thick ice point cloud subset. The surface normal vector is calculated based on the thick ice point cloud subset to obtain the normal vector set; The set of normal vectors is compared with the preset reference contour normal vectors, local extreme points are filtered out, and aggregated to obtain the ice layer protrusion region; Density analysis was performed on the ice protrusion region to obtain the spatial distribution density matrix; The volume integral of the ice protrusion region is performed based on the spatial distribution density matrix to obtain the first ice layer thickness as the ice layer thickness data.

[0041] The spatial distribution density matrix records the number of point clouds within each 3D grid and is used to filter grids containing valid ice layer point clouds. The local heights of each grid are arranged separately according to their grid positions to form a local height matrix. The spatial distribution density matrix is ​​used to determine valid grids, and the local height matrix is ​​used for subsequent volume integration calculations. The local height of the ice layer in each valid grid is determined separately based on the distance along the local normal direction between points on the outer surface of the ice layer and the ice-free reference surface.

[0042] In one possible implementation, if the icing scale span is less than or equal to a preset scale span threshold, then the frost thickness is inverted on the ice point cloud data to obtain a second ice layer thickness as ice layer thickness data, including: If the ice scale span is less than or equal to the scale span threshold, feature extraction is performed on the ice layer point cloud data to obtain the point cloud roughness matrix; Gradient frequency domain transformation is performed on the point cloud roughness matrix to obtain the texture frequency matrix; Element-wise analysis of the texture frequency matrix yields the frost coverage rate; The refractive offset is obtained by refraction derivation based on the frost coverage rate. The refraction offset is used to perform refraction inversion calculations to obtain the second ice layer thickness as ice layer thickness data.

[0043] It should be noted that the preset scale span threshold is fixed at 25 mm. This threshold was determined based on classification and verification of 300 sets of blade icing area samples labeled with thin frost and thick ice categories. Using a span of 25 mm as the boundary, the best distinction between thin frost and thick ice samples can be achieved, with a verification accuracy of 93%. The average ice thickness corresponding to this threshold is approximately 20 mm, which is basically consistent with the icing level that requires de-icing operations in wind farm operation and maintenance, and is the critical value for distinguishing between thick ice accumulation and thin frost adhesion. The icing scale span is a single value calculated for the entire icing area divided by the initial icing boundary. The icing scale span is compared with the scale span threshold. If the icing scale span is greater than the threshold, it is determined that there is severe thick ice accumulation in the area, and the thick ice thickness calculation branch is initiated. The thick ice point cloud subset is clipped based on the spatial boundary of the geometric convex set. All point clouds within the corresponding range are clipped from the original ice point cloud data to obtain the thick ice point cloud subset. Subset point clouds are specifically designed for fine-grained thickness calculations in thick ice regions. They can eliminate interference from surrounding thin ice regions, and the point cloud size is approximately 20% to 50% of the original point cloud, which can improve the calculation efficiency and accuracy of thick ice thickness.

[0044] It is worth noting that the surface normal vector of each data point within the thick ice point cloud subset is calculated using the same neighborhood covariance method as in step S13. That is, the covariance matrix is ​​constructed by taking the 20 nearest neighbors, and the eigenvector corresponding to the smallest eigenvalue is the normal vector of that point. The direction is uniformly pointed towards the outer side of the blade. The reference contour normal vector of the standard 3D model of the blade in ice-free condition is pre-stored, with corresponding reference normal vector data for each spatial location. The normal vector of the thick ice point cloud is compared point-by-point with the reference contour normal vector, and the spatial angle between the two sets of normal vectors is calculated using the vector dot product angle formula. If the angle is greater than 45 degrees, the point is determined to be a local extremum. This angle threshold is based on the statistical distribution of normal vectors from 50 sets of measured ice nodule samples. The angle between the normal vector of the protruding area and the reference normal is generally greater than 40 degrees; a 45-degree threshold ensures a detection rate of over 95% for ice protrusions. Spatially adjacent local extremum points are aggregated to form continuous connected regions, which are the ice protrusion regions, corresponding to thick ice core structures such as ice nodules and ice ridges. The neighborhood distance threshold is set to 3 mm. This neighborhood distance threshold is based on the statistically typical size of ice nodules and ice ridges on the blade surface. The distance between adjacent independent protrusions is usually greater than 5 mm; the 3 mm threshold can merge extremum points of the same ice nodule while avoiding the erroneous merging of different ice ridge structures.

[0045] It is worth further explaining that the calculation of the first ice layer thickness is based on volume integration. The ice protrusion region is divided into a three-dimensional grid with a side length of 5 mm. The number of point clouds in each grid is counted, and the grids are arranged according to their spatial positions to obtain a spatial distribution density matrix. Each element in the spatial distribution density matrix represents the number of point clouds in the corresponding grid. Based on the spatial distribution density matrix, the effective grids with continuous point cloud distribution are determined. For each effective grid, the local normal direction of the blade's standard three-dimensional model at the corresponding position in that grid is used as the thickness direction. The distance between the points on the outer surface of the ice layer in the grid and the ice-free reference surface is calculated, and the average of these distances is taken as the local ice layer height of that grid. The local ice layer height of each effective grid is multiplied by the projected area of ​​the corresponding grid on the ice-free reference surface to obtain the local ice layer volume corresponding to each effective grid. All local ice layer volumes are accumulated to obtain the estimated total volume of the ice protrusion region. The estimated total volume is divided by the projected area of ​​all effective grids to obtain the initial average thickness. This initial average thickness is then corrected using a pre-calibrated dimensionless thickness calibration coefficient to obtain the first ice layer thickness, measured in millimeters. The thickness calibration coefficient is calibrated based on the correspondence between the point cloud thickness of a standard ice layer sample and the measured thickness.

[0046] Subsequently, if the ice scale span is less than or equal to the scale span threshold, the corresponding area is determined to be in a state of slight frost adhesion, and the frost thickness calculation branch is initiated. In this type of scenario, the geometric protrusions of the ice layer are not obvious, and directly calculating the thickness using point cloud data would result in a large error. Therefore, a method combining topographic roughness and optical properties is used to invert the thickness. The extraction of the point cloud roughness matrix is ​​based on the calculation of local geometric undulations of the point cloud. First, the ice layer point cloud data is subjected to two-dimensional rasterized projection, with the raster size corresponding to the pixels of the original image. The elevation standard deviation of all points within each raster is calculated, and this standard deviation represents the surface roughness. Then, the standard deviation is linearly mapped to a gray level of 0-255 to form a two-dimensional matrix of the same size as the original image, which is the point cloud roughness matrix. The higher the gray value, the more obvious the surface geometric undulations, and the denser the corresponding frost crystal structure. Unlike image-based roughness matrices, the point cloud roughness matrix is ​​derived from the real geometric topography of the three-dimensional point cloud, is not affected by visual interference such as lighting or stains, and can more accurately reflect the surface micro-features in the state of thin frost. The spatial resolution of the point cloud roughness matrix matches the density of the original point cloud, which is 1 pixel per millimeter, providing basic data for subsequent frost coverage calculation.

[0047] It is worth noting that, for gradient frequency domain transformation based on the point cloud roughness matrix, the gradient magnitude and direction of each pixel must first be calculated using the Sobel operator to obtain a gradient magnitude map. Then, a two-dimensional fast Fourier transform is performed on the gradient magnitude map to obtain a spectrum map. The coefficients in the low-frequency region at the center of the spectrum are set to zero to achieve high-pass filtering, preserving high-frequency components. Next, the high-frequency energy within each 8×8 sub-block of the filtered spectrum map is statistically analyzed to construct a dimension-reduced texture frequency matrix, where the matrix elements are the logarithmic values ​​of the high-frequency energy in each sub-block. Then, element-wise analysis is performed on the texture frequency matrix to obtain a pre-formed reference texture frequency matrix for the same leaf region under frost-free conditions. The current texture frequency matrix is ​​compared element-wise with the reference texture frequency matrix, and the proportion of elements with current high-frequency energy higher than the corresponding reference high-frequency energy is counted out of the total number of elements to obtain the frost coverage rate. Before calculating the refraction offset, the current thin frost image and the frost-free reference image are first geometrically registered using the fixed texture features of the leaves outside the icy area to eliminate the overall coordinate deviation caused by camera pose changes and leaf displacement. Then, the corresponding texture feature points inside the icy area are matched, and the average value of the local residual coordinate offsets after registration is used as the refraction offset. This value reflects the image pixel displacement caused by refraction after light passes through the frost layer. Subsequently, a pre-calibrated "frost coverage - refraction offset - frost thickness" lookup table is used, with frost coverage and refraction offset as inputs, to output the second ice layer thickness in millimeters via bilinear interpolation. The lookup table is used in a controlled-temperature environment of -5℃ to perform transmission calibration on frost samples in the 0 to 20 mm range with a step size of 0.5 mm. Standard frost samples with a thickness of 0 to 20 mm are prepared with a step size of 0.5 mm. The true value of the frost thickness of each sample is measured using a laser displacement sensor. At the same time, the vision module acquires the icing image of each sample and extracts the frost coverage and refractive offset. The three sets of data are matched one-to-one to construct a two-dimensional lookup table. For example, in a thin frost area, the frost coverage is 65% and the measured refractive offset is 2.8 pixels. After interpolation using the lookup table, the second ice layer thickness is found to be 0.8 mm.

[0048] In step S15, the ice thickness data is classified using multidimensional features to obtain the degree of icing. Based on the degree of icing and pre-acquired historical icing data, a difference calculation is performed to obtain the icing evolution rate, including: The thickness difference matrix is ​​obtained by calculating the pixel-by-pixel difference of the ice layer thickness data; Multiple features are extracted from the thickness difference matrix to obtain a fused feature vector; The fused feature vectors are mapped to a high-dimensional hyperplane to obtain the decision boundary distance; The degree of icing is obtained by classifying the data based on the decision boundary distance. Based on the degree of icing, obtain historical icing data for the corresponding collection period; The icing evolution rate is obtained by extrapolating the icing degree and the historical icing data.

[0049] It should be noted that the ice thickness data in this step is uniformly converted according to the same spatial coordinates, normal direction, and thickness measurement reference as the blade thickness matrix when the blade is ice-free, before participating in the pixel-by-pixel difference calculation. Specifically, the ice thickness of each pixel position obtained in step S14 is superimposed on the ice-free blade thickness reference corresponding to the same pixel position to obtain the total thickness value of the corresponding pixel position under the current icing state, and the total thickness value is used as the ice thickness data participating in the difference calculation in this step. Subtracting the ice-free blade thickness value at the same pixel position from the total thickness value, the resulting difference is the increase in ice thickness at that pixel position relative to the ice-free state. Therefore, the pixel-by-pixel difference calculation does not directly subtract the ice thickness from the blade body thickness, but rather obtains the ice thickness increase by subtracting the current total thickness under the same measurement reference from the ice-free reference thickness.

[0050] It should be noted that the thickness difference matrix is ​​obtained by subtracting the ice thickness data from the blade thickness matrix when the blade is ice-free, pixel by pixel. The blade thickness matrix when the blade is ice-free is derived from the standard 3D model of the blade introduced in step S13. This 3D model establishes a coordinate system with the center of the flange at the blade root as the origin, the spanwise direction as the X-axis, the chord direction as the Y-axis, and the thickness direction as the Z-axis. The surface accuracy of the model is ±1 mm. When generating the design thickness matrix, the same spatial resolution and pixel coordinate system as the ice thickness data matrix are used. The standard 3D model of the blade is projected onto the 2D imaging plane along the thickness direction, and the normal distance between the upper and lower surfaces of the blade is taken as the design thickness value at each pixel location. The ice thickness data is stored in the form of a pixel-level spatial matrix with the dimension of millimeters. Each element in the matrix corresponds to the ice thickness value at the corresponding pixel location. The ice thickness data can intuitively reflect the spatial distribution gradient and thickness abrupt boundary of the ice thickness, containing rich information such as thickness fluctuations and edge abrupt changes, providing a basic data source for subsequent multi-feature extraction.

[0051] Subsequently, multi-feature extraction includes two types of features: spatial distribution features and edge contour gradients. Spatial distribution features are obtained by statistically analyzing the gray-level distribution of the thickness matrix, extracting two indicators: fluctuation amplitude and uniformity, to quantify the thickness fluctuation amplitude and distribution uniformity within the matrix, respectively. Edge contour gradients are obtained by detecting the boundary positions of abrupt thickness changes, outputting feature information such as gradient magnitude, length, and density. Both types of features are normalized to the 0-1 range to eliminate dimensional differences before linear concatenation, generating a one-dimensional fused feature vector with a total length of 128 dimensions. The first 64 dimensions correspond to multi-scale spatial distribution features, and the last 64 dimensions correspond to edge gradient features. The feature extraction window size and computation step size are calibrated using samples to control computational load while ensuring feature representation capabilities. The fused feature vector can comprehensively represent the spatial distribution pattern and boundary characteristics of ice layer thickness, providing multi-dimensional input for the classification model, avoiding classification bias caused by single features, and improving the stability of subsequent icing degree determination.

[0052] It is worth noting that the high-dimensional hyperplane mapping is achieved through a pre-trained classifier, which maps the fused feature vectors to a high-dimensional feature space, constructing an optimal classification hyperplane in the high-dimensional space that maximizes the margin between the two classes of samples. The perpendicular geometric distance from the feature vector mapping point to the optimal classification hyperplane is calculated; this distance is the decision boundary distance. The classifier is trained on 1500 sets of labeled ice thickness samples, covering all levels of icing from thin frost to thick ice, and undergoes multiple rounds of cross-validation to ensure generalization ability. The decision boundary distance is a scalar value; a positive value indicates that the feature vector falls within the thick ice half-space, and a negative value indicates that it falls within the thin ice half-space. The larger the absolute value of the decision boundary distance, the higher the confidence of the classification decision. Specifically, the mapping from the decision boundary distance to the degree of icing is as follows: when... At that time, the confidence level for thin ice was extremely high, based on the numerical value of the degree of freezing. Mapped to levels 1-5; when At that time, the confidence level for thick ice was extremely high, based on the numerical value of the degree of freezing. Mapped to levels 6-10. When At this point, insufficient classification confidence and samples close to the classification boundary trigger the deep physical feature extraction process. It's worth further explaining that the classification process first compares the absolute value of the decision boundary distance with a preset classification threshold of 2. This threshold is based on 5-fold cross-validation using 1500 labeled samples, and the decision boundary distance value of 2 corresponding to a 95% classification confidence level is used as the threshold, ensuring classification accuracy under normal icing conditions. If the absolute value of the decision boundary distance... If the classification result has sufficient confidence, the corresponding icing level can be directly output based on the decision boundary distance d; if the absolute value of the decision boundary distance is... This indicates a unique icing condition, triggering a deep physical feature extraction process to further improve classification accuracy. Deep features include two types of data: ice crystal density distribution and structural stress distribution. Ice crystal density distribution is obtained by collecting and calculating data using ultrasonic detection units deployed inside the blade's leading edge. These units are arranged at 2-meter intervals along the blade's span, with a sampling frequency of 1 Hz. During the acquisition process, ultrasonic signals are emitted into the ice layer, and the propagation time of the signals through the ice is recorded. The propagation speed of the sound waves within the ice layer is calculated by combining this with the ice layer thickness at the corresponding location. The ice crystal density value at each measuring point is then calculated by comparing it with a pre-calibrated relationship between sound velocity and density. The density value ranges from 700 to 900 kg / m³, with low density values ​​corresponding to frost-like ice and high density values ​​corresponding to transparent ice. The density calculation error for a single measuring point is controlled within ±30 kg / m³. Structural stress distribution is obtained by collecting deformation data using strain detection units attached to the blade's leading edge surface, combined with blade structural parameters and the ice layer load distribution pattern. These strain detection units are evenly spaced along the chord direction, with a sampling frequency of 10 Hz. The calculation process uses measured deformation as boundary conditions, deriving structural stress values ​​induced by ice adhesion region by region, and outputting two types of distribution parameters: shear stress and normal stress, all in megapascals (MPA). The relative error of the stress calculation results does not exceed 8%, accurately reflecting the changes in blade structural loads caused by ice adhesion. A comprehensive classification is performed combining these two types of deep physical characteristics. Ice crystal density, shear stress, and normal stress are normalized to their maximum and minimum values, and a weighted comprehensive score is calculated with ice crystal density weighting 0.5, shear stress weighting 0.3, and normal stress weighting 0.2, ranging from 0 to 1. The comprehensive score is multiplied by 9, then added by 1, rounded, and linearly mapped to a level from 1 to 10 to obtain the degree of icing. Higher levels indicate more severe icing; for example, severe icing corresponds to level 8.

[0053] It should be noted that historical icing data is stored in a local time-series database, consisting of icing severity values ​​from past sampling periods bound to corresponding timestamps. Each data point also includes a corresponding blade region number. The data acquisition period is fixed at 5 minutes by default. This period is set based on the average evolution rate of icing under low-temperature conditions, balancing real-time monitoring and data processing efficiency. When the ambient temperature is below -10 degrees Celsius, it can be automatically shortened to 2 minutes. When acquiring historical icing data, the current icing severity acquisition timestamp is used as the baseline, and continuous sampling data from the past 2 hours are extracted to construct an icing severity time series. The data in the series is arranged in ascending order of time, and each data point is bound to a precise second-level timestamp. Validity verification is performed simultaneously during the data extraction process. The preset anomaly judgment difference is 8 levels, which is set based on the icing severity level range and the sampling period. Under the default sampling period of 5 minutes, the icing severity level range is 1 to 10. An icing severity difference exceeding 8 levels between two consecutive sampling periods is physically impossible and is judged as sensor failure or abnormal data identification, and is replaced using linear interpolation. The anomaly difference setting is independent of the rate threshold. The rate threshold is preset to one level per minute, corresponding to a difference threshold of five levels within a 5-minute sampling period. The anomaly judgment difference of eight levels is much larger than the rate judgment difference of five levels, ensuring that rapid freezing signals are not deleted during the validity verification stage, thus achieving parallel processing of anomaly data filtering and rapid freezing early warning. Field operation statistics show a data effectiveness rate of no less than 99.5%, and the complete time series can present the evolution trend of freezing degree over time.

[0054] The evolution calculation first extracts the icing degree values ​​from two adjacent collection periods in the time series. The value of the later period is subtracted from the value of the earlier period to obtain the icing degree difference. A positive difference indicates increased icing, while a negative difference indicates melting. The scheme focuses on risk scenarios of increased icing. A preset incremental threshold of level 2 is used. This threshold is statistically calibrated based on 200 sets of field monitoring data under normal icing conditions from -5℃ to 0℃. The average icing degree increment within 5 minutes is level 0.08, with a standard deviation σ = level 0.08. The threshold level 2 is obtained by adding 25 times the standard deviation to the mean (i.e., 0.08 + 25 × 0.08 = 2.08) and rounding down, which can effectively identify abnormally rapid icing processes. If the difference in icing intensity is less than or equal to the increment threshold, it is determined to be normal slow icing, and the evolution rate is calculated as the quotient of the difference and the time span. If the difference is greater than the increment threshold, the time span corresponding to two cycles is further extracted, and the quotient is calculated using the difference in icing intensity as the dividend and the time span as the divisor to obtain the icing evolution rate, measured in levels per minute. For example, if the icing intensity increases from level 3 to level 8 within 5 consecutive minutes, the difference is 5 levels, corresponding to an evolution rate of 1 level per minute. This rate can quantify the speed at which icing intensifies.

[0055] In step S16, the evolution duration is extracted based on the icing evolution rate, and the evolution type is matched based on the evolution duration to obtain the final icing state, including: Based on the freezing evolution rate and the preset rate threshold, time series analysis is performed to obtain the evolution duration. The physical evolution type is obtained by matching the evolution type based on the duration of the evolution. The final freezing state is obtained by mapping the state according to the physical evolution type.

[0056] It should be noted that the preset rate threshold is set to level 1 per minute. This threshold is determined based on statistical data of wind turbine icing accidents and meteorological observation records from the past three years. In scenarios where rapid freezing leads to blade damage, the icing evolution rate is generally greater than level 0.8 per minute. A safety factor of 1.25 is used, rounded up to obtain level 1.0 per minute as the final threshold, which is the critical rate distinguishing between conventional icing and rapid freezing. The icing evolution rate is compared with the rate threshold. If the icing evolution rate is greater than the rate threshold and the condition is met for two consecutive sampling periods, it is determined that a rapid freezing state has been entered. At the same time, the process of extracting the evolution duration is triggered to avoid misjudgment caused by fluctuations in single-period data. The evolution duration is calculated through time series backtracking. Starting from the timestamp of the current moment, the time series is traced backward until the first time node that is continuously lower than the rate threshold is found. This node is the starting point of the high evolution rate mutation. The time difference between the current moment and the mutation point is calculated. This difference is the evolution duration, measured in minutes. For example, a high evolution rate started at 9:10 AM, and the current time is 9:30 AM, corresponding to an evolution duration of 20 minutes.

[0057] It is worth noting that the evolution type matching is based on preset classification rules. The input consists of two features: icing evolution rate and evolution duration, and the output is the corresponding physical evolution type. The classification rules are based on statistics from 300 historical icing process samples, covering various icing evolution processes under different temperatures, humidity, and wind speeds, including three categories: rapid freezing, slow icing, and icing melting. The classification accuracy can reach over 92%. Among the two features, the icing evolution rate accounts for approximately 60% of the weight, and the evolution duration accounts for approximately 40%. The classification rules clearly state that when the icing evolution rate is higher than the rate threshold and the evolution duration exceeds 10 minutes, the corresponding process is determined to be a rapid freezing physical evolution; if the icing evolution rate meets the threshold but the duration is less than 10 minutes, it is determined to be a short-term, small fluctuation, not a stable physical evolution process. Rapid freezing physical evolution corresponds to the rapid impact and freezing of supercooled water droplets or the rapid growth of ice layers caused by a sudden drop in temperature. This type of process has a significant impact on wind turbine load and is a high-risk condition that is the focus of monitoring.

[0058] It is worth further explaining that the state mapping first converts the physical evolution type into a corresponding icing state identifier, using integer encoding to distinguish different states. For example, 1 represents normal light icing, 2 represents slowly thickening icing, 3 represents moderate-risk icing, and 4 represents rapid freezing with severe risk. Each code corresponds to a specific operation and maintenance handling strategy. After the state identifier is generated, the final icing state of the blades is comprehensively determined by combining the design load-bearing limit of the wind turbine blades with the current icing degree. When the rapid freezing type is combined with an icing degree of level 8 or above, it is determined to be a severe-risk icing state, requiring immediate triggering of de-icing operations and load reduction protection. The final icing state includes five core information categories: spatial coordinates of the icing area, average ice thickness, maximum ice thickness, evolution type, and risk level. The data is transmitted to the wind turbine main control system through an industrial standard protocol, serving as the basis for decision-making in de-icing strategy formulation and operational parameter adjustment. A pilot verification at a mountain wind farm showed that the accuracy rate of this rapid freezing early warning can reach over 90%, effectively reducing the risk of equipment damage caused by icing and improving the safety and reliability of wind turbine operation in winter.

[0059] In summary, this invention discloses a method for identifying icing regions on wind turbine blades, including acquiring an original image of the wind turbine blade; performing illumination equalization correction on the original image to obtain a corrected image; extracting texture features from the corrected image to obtain a roughness matrix; performing region segmentation based on the roughness matrix to obtain an initial icing boundary; extracting internal point clouds based on the initial icing boundary to obtain ice layer point cloud data; extracting geometric protrusions based on the ice layer point cloud data to obtain the icing scale span; if the icing scale span is greater than a preset scale span threshold, then the ice layer point cloud is segmented. The cloud data is used to calculate the thickness of the first ice layer as the ice layer thickness data. If the icing scale span is less than or equal to a preset scale span threshold, the frost thickness is inverted from the ice layer point cloud data to obtain the second ice layer thickness as the ice layer thickness data. The ice layer thickness data is then classified using multi-dimensional features to obtain the degree of icing. The difference between the degree of icing and pre-acquired historical icing data is calculated to obtain the icing evolution rate. The evolution duration is extracted based on the icing evolution rate, and the evolution type is matched based on the evolution duration to obtain the final icing state. This invention achieves accurate identification of icing areas on wind turbine blades under complex operating conditions through image illumination correction, blade surface texture analysis, and fusion of three-dimensional morphological features.

[0060] The second embodiment of the present invention provides a wind turbine blade icing area identification system, including: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the steps of the method described above.

[0061] It should be noted that the wind turbine blade icing area identification system provided in this embodiment of the invention is used to execute all the process steps of the wind turbine blade icing area identification method in the above embodiment. The working principle and beneficial effects of the two are one-to-one, so they will not be described again.

[0062] It should be noted that the device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Furthermore, in the accompanying drawings of the device embodiments provided by this invention, the connection relationships between modules indicate that they have communication connections, which can be specifically implemented as one or more communication buses or signal lines. Those skilled in the art can understand and implement this without any creative effort.

[0063] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A method for identifying icing areas on wind turbine blades, characterized in that, include: Acquire the original image of the wind turbine blades, and perform illumination equalization correction on the original image to obtain the corrected image; Texture features are extracted from the corrected image to obtain a roughness matrix. Region segmentation is then performed based on the roughness matrix to obtain the initial icing boundary. Based on the initial icing boundary, internal point cloud is extracted to obtain ice layer point cloud data. Based on the ice layer point cloud data, geometric protrusions are extracted to obtain the icing scale span. If the ice scale span is greater than a preset scale span threshold, the ice layer point cloud data is used to calculate the ice thickness to obtain the first ice layer thickness as the ice layer thickness data. If the icing scale span is less than or equal to a preset scale span threshold, then the frost thickness is inverted on the ice point cloud data to obtain the second ice layer thickness as the ice layer thickness data. The ice thickness data is classified using multidimensional features to obtain the degree of icing. The difference between the degree of icing and the pre-acquired historical icing data is calculated to obtain the icing evolution rate. The duration of the icing evolution is extracted based on the icing evolution rate, and the evolution type is matched based on the duration of the icing evolution to obtain the final icing state.

2. The method for identifying icing areas on wind turbine blades according to claim 1, characterized in that, The process of acquiring the original image of the wind turbine blades and performing illumination equalization correction on the original image to obtain a corrected image includes: Acquire the original image of the wind turbine blade, extract the original pixel matrix of the original image, and generate a grayscale distribution histogram; Illumination distance analysis is performed based on the grayscale distribution histogram to obtain the mapping function; The target pixel matrix is ​​obtained by non-linearly stretching the original pixel matrix using the mapping function. The image is reconstructed based on the target pixel matrix to obtain the corrected image.

3. The method for identifying icing areas on wind turbine blades according to claim 1, characterized in that, The process of extracting texture features from the corrected image to obtain a roughness matrix, and then performing region segmentation based on the roughness matrix to obtain the initial icing boundary, includes: Pixel-to-pixel frequency statistics are performed on the corrected image to obtain the gray-level co-occurrence matrix; Feature extraction is performed based on the gray-level co-occurrence matrix to obtain gray-level image features, and spatial mapping is performed on the gray-level image features to obtain local roughness; The roughness matrix is ​​obtained by spatially arranging the local roughness. The roughness matrix is ​​used to predict icing, generating an icing probability matrix; The initial icing boundary is obtained by performing boundary segmentation based on the icing probability matrix.

4. The method for identifying icing areas on wind turbine blades according to claim 1, characterized in that, The step of extracting internal point clouds based on the initial icing boundary to obtain ice layer point cloud data, and extracting geometric protrusions based on the ice layer point cloud data to obtain the icing scale span, includes: The area within the initial icing boundary is extracted to obtain the ice layer region. The coordinates of the ice layer region are then transformed to obtain ice layer point cloud data. Curvature filtering is performed on the ice point cloud data to obtain a set of geometric protrusions; Calculate the vertical distance between the set of geometric protrusions and the preset reference plane to obtain the height dataset; The maximum height difference is calculated based on the height dataset to obtain the icing scale span.

5. The method for identifying icing areas on wind turbine blades according to claim 1, characterized in that, If the icing scale span is greater than a preset scale span threshold, then the ice layer point cloud data is used to calculate the ice thickness to obtain a first ice layer thickness as ice layer thickness data, including: If the ice scale span is greater than a preset scale span threshold, the ice point cloud data is cropped to obtain a thick ice point cloud subset. The surface normal vector is calculated based on the thick ice point cloud subset to obtain the normal vector set; The set of normal vectors is compared with the preset reference contour normal vectors, local extreme points are filtered out, and aggregated to obtain the ice layer protrusion region; Density analysis was performed on the ice protrusion region to obtain the spatial distribution density matrix; The volume integral of the ice protrusion region is performed based on the spatial distribution density matrix to obtain the first ice layer thickness as the ice layer thickness data.

6. The method for identifying icing areas on wind turbine blades according to claim 1, characterized in that, If the icing scale span is less than or equal to a preset scale span threshold, then the frost thickness is inverted on the ice point cloud data to obtain a second ice layer thickness as ice layer thickness data, including: If the ice scale span is less than or equal to the scale span threshold, feature extraction is performed on the ice layer point cloud data to obtain the point cloud roughness matrix; Gradient frequency domain transformation is performed on the point cloud roughness matrix to obtain the texture frequency matrix; The frost coverage rate is obtained by performing element-wise statistics on the texture frequency matrix. The refractive offset is obtained by refraction derivation based on the frost coverage rate. The refraction offset is used to perform refraction inversion calculations to obtain the second ice layer thickness as ice layer thickness data.

7. The method for identifying icing areas on wind turbine blades according to claim 1, characterized in that, The process of classifying the ice thickness data using multidimensional features to obtain the degree of icing, and calculating the difference between the degree of icing and pre-acquired historical icing data to obtain the icing evolution rate, includes: Multiple features are extracted from the ice thickness data to obtain a fused feature vector; The fused feature vectors are mapped to a high-dimensional hyperplane to obtain the decision boundary distance; The degree of icing is obtained by classifying the data based on the decision boundary distance. Based on the degree of icing, obtain historical icing data for the corresponding collection period; The icing evolution rate is obtained by extrapolating the icing degree and the historical icing data.

8. The method for identifying icing areas on wind turbine blades according to claim 1, characterized in that, The step of extracting the evolution duration based on the freezing evolution rate, performing evolution type matching based on the evolution duration, and obtaining the final freezing state includes: Based on the freezing evolution rate and the preset rate threshold, time series analysis is performed to obtain the evolution duration. The physical evolution type is obtained by matching the evolution type based on the duration of the evolution. The final freezing state is obtained by mapping the state according to the physical evolution type.

9. A wind turbine blade icing area identification system, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the method according to any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor, it implements the method described in any one of claims 1 to 8.

Images